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ABSTRACT 


Detecting  short  transients  is  a  signal  processing  application  that  has  a  wide  range 
of  military  uses.  To  be  specific  in  Undersea  Warfare,  sensitive  signal  detection  schemes 
can  increase  the  effective  range  of  active  and  passive  sonar  operations.  Current  research 
is  being  done  to  improve  the  capability  of  detecting  short  signals  buried  within 
background  noise,  particularly  in  Littoral  waters.  Starting  with  a  colored  noise  model, 
this  thesis  will  introduce  two  denoising  methods  based  on  multiresolution  analysis  and 
compare  the  results  to  current  transient  detection  techniques.  The  goal  of  this  thesis  is 
not  necessarily  to  replace  current  detection  schemes,  but  rather  to  enhance  them  and 
thereby  making  the  procedure  more  robust. 
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I. 


INTRODUCTION 


A.  PURPOSE 

Since  the  end  of  the  Cold  War,  there  has  been  a  paradigm  shift  in  detecting, 
tracking,  classifying  and  ultimately  neutralizing  subsurface  threats.  In  its  most  general 
form,  the  US  Navy  now  sees  a  more  probable  warfare  scenario  in  coastal  waters  as 
opposed  to  the  Blue  Water  confrontations  occurring  earlier  last  century.  The  major  threat 
of  the  once  large  Soviet  Nuclear  Submarine  Force  has  now  turned  into  a  small-country 
coastal  diesel  submarine  threat.  Indeed,  the  very  name  for  this  warfare  area  changing 
from  Anti-submarine  Warfare  to  Undersea  Warfare  implies  a  change  in  our  subsurface 
warfare  strategy. 

With  the  shift  in  subsurface  threats  come  different  challenges  in  detection  and 
classification.  No  longer  can  acoustic  detection  systems  exploit  deep  sound  channels  on 
a  relatively  noisy  nuclear  submarine  threat.  In  this  new  era  where  the  threat  of  the  once 
large  nuclear  fleet  from  the  Former  Soviet  Union  is  limited  at  best,  a  new  foe  has 
emerged  in  coastal  waters.  Not  only  do  Naval  USW  systems  have  to  contend  with 
increased  challenges  such  as  bottom  bounce  and  added  noise  due  to  a  higher  density  of 
neutral  shipping,  but  the  platforms  of  choice  are  the  quiet  diesel  submarines  while  on 
battery  power.  In  addition  to  this  diesel  threat  in  a  Littoral  Warfare  environment,  mines 
have  become  more  of  a  threat.  Because  of  their  high  destruction-to-cost  ratio,  laying 
mines  to  thwart  US  global  interests  is  a  cost-effective  way  to  inflict  maximum  damage 
with  little  cost.  Indeed,  the  US  Navy  has  suffered  more  battle  damage  to  their  ships  due 
to  mines  than  any  other  subsurface  threat  since  the  Korean  War.  Mines  are  cheap,  easy  to 
deploy  and  hard  to  detect.  To  many  countries,  mines  are  a  viable  substitute  for  an 
otherwise  weak  coastal  defense. 

The  nature  of  Littoral  Warfare  possesses  unique  challenges  to  detecting  and 
classifying  subsurface  threats.  Since  the  Walker  Treason  in  the  middle  1980’s,  enemy 
acoustic  signatures  have  been  suppressed  to  an  all-time  low.  One  way  to  regain  a 
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considerable  edge  in  Undersea  Warfare  is  to  investigate  more  robust  methods  of  detecting 
enemy  acoustic  signatures  within  heavy  background  noise. 


B.  GOALS 

This  thesis  is  part  of  a  larger  project  that  addresses  the  issue  of  improving  the 
tools  available  to  USW  systems  operators  of  surface  and  subsurface  platforms. 
Currently,  USW  systems  operators  use  variations  of  a  spectrogram  as  their  main  visual 
display  to  monitor  passive  acoustic  signatures.  This  type  of  analysis  is  good  for  signals 
that  are  long  in  duration.  Existing  displays  show  a  long-term  history  of  acoustic  data, 
which  is  good  for  determining  underwater  activity  over  a  long  period  of  time.  Using  the 
same  techniques,  short-duration  signals  (or  transients),  however,  may  go  unnoticed.  For 
example,  if  an  operator  is  either  fatigued  or  focusing  on  long-term  trends,  he  might 
overlook  a  very  slight  indication  that  a  transient  signal  a  fraction  of  a  second  in  duration 
had  just  occurred. 

The  thesis  goal  is  to  create  more  robust  visual  and  audible  features  at  the 
operator’s  disposal  to  enhance  the  detection  of  transient  signatures.  Specifically,  the 
thesis  will  develop  and  investigate  three  signal  processing  procedures  based  on  Wavelet 
and  Multiresolution  analyses,  which  will  attempt  to  detect  transient  signals  currently 
overlooked  by  common  signal  processing  techniques. 

This  thesis  will  also  investigate  an  automation  procedure,  which  will  alarm  the 
operator  of  such  transient  signals. 


C.  METHODOLOGY 

It  is  the  main  thrust  of  this  thesis  to  detect  transient  signals  in  a  noisy  acoustic 
environment.  One  way  to  investigate  the  improvement  of  detecting  desired  acoustic 
signatures  imbedded  within  a  noisy  environment  is  by  studying  the  statistics  and  nature 
of  the  signals.  Although  littoral  warfare  might  be  difficult  in  terms  of  background 
shipping  and  silent  threats,  one  positive  aspect  is  that  acoustic  data  is  easily  accessible 
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because  the  surveillance  area  is  small  and  there  exist  well-known  shipping  lanes  and 
underwater  topographies. 

Current  methods  of  transient  detection  use  time-frequency  analysis  and  stochastic 
modeling.  Although  these  methods  have  been  proven  quite  useful,  they  have  their 
limitations.  This  thesis  will  investigate  the  root  causes  of  these  weaknesses  and  explore 
some  ideas  to  circumvent  them.  Ideas  such  as  Filter  Banks  using  Wavelet  Transforms 
and  Multiresolution  Filtering  just  might  be  the  missing  tools  to  enhance  transient 
detection. 

Three  methods  have  been  developed  to  detect  transient  signals  in  noise.  One  of 
these  methods  is  based  on  Wavelet  Analysis  and  the  other  two  based  on  an  ARMA  model 
and  multiresolution  processing.  A  standard  data  set  is  used  to  compare  the  new  signal 
processing  methods  against  well-established  methods.  From  this  data,  conclusions  will 
be  drawn  and  recommendations  will  be  made  for  further  improvement  of  acoustic 
transient  detection  siystems. 


D.  BENEFITS 

The  confidence  of  a  decision  is  directly  proportional  to  the  amount  of  concurring 
evidence.  Consider  the  method  of  detecting  submarines  using  the  classic  method  of 
Target  Motion  Analysis  (TMA).  Since  TMA  uses  nothing  but  passive  techniques, 
determining  a  target’s  bearing,  range,  course  and  speed  can  be  challenging.  What 
strengthens  this  method  of  detecting  underwater  weapons  systems  is  the  redundancy  of 
evidence.  Information  that  does  not  show  up  on  a  Dead  Reckoning  Trace  (DRT)  plot 
might  be  found  on  the  time-frequency  table  or  on  a  Moboard  plot.  Given  a  well-behaved 
target,  the  success  of  TMA  relies  on  the  fact  that  dynamic  changes  will  be  detected  in 
different  domains.  For  example,  a  shift  in  frequency  may  be  the  first  detection  of  a 
change  in  course  or  speed,  but  probably  won’t  be  the  first  method  to  estimate  a  target’s 
range  or  course.  The  different  techniques  in  TMA  are  independent  in  the  sense  that  they 
do  not  rely  on  similar  assumptions.  These  methods  each  provide  a  piece  of  a  puzzle, 
which  when  put  together,  can  accurately  determine  particular  warfare  scenario. 
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Like  the  set  of  tools  used  in  TMA,  an  acoustic  transient  detection  system  needs  to 
be  developed  which  takes  the  same  fact-layering  approach  in  order  to  make  an  informed 
decision.  This  thesis  does  not  attempt  to  improve  on  well-established  signal  processing 
techniques  such  as  time-frequency  analysis.  Rather,  the  goal  of  this  thesis  is  to  create 
additional  layers  (based  on  fundamentally  different  assumptions)  and  add  to  the  fact¬ 
finding  set  in  order  to  enhance  the  robustness  of  detecting  transient  signals. 

The  benefits  of  this  thesis  will  be  the  addition  of  new  tools,  which  will  enhance  a 
transient  detection  system.  One  of  the  new  signal  processing  techniques  may  be  the 
necessary  tool  needed  to  detect  some  transient  signals. 

Another  benefit  to  the  undersea  warfighting  capability  of  a  ship  is  to  discuss 
methods  to  enhance  the  automation  of  detecting  valid  transient  signals. 


4 


II.  CURRENT  TRANSIENT  SIGNAL  DETECTION  TOOLS 

Although  current  methods  of  detection  can  successfully  identify  transient  signals, 
there  is  a  threshold  floor  or  frequency  bandwidth  to  which  signals  become 
indistinguishable  from  background  noise.  Of  those  cases  where  current  detection 
schemes  miss  a  signal,  we  will  show  that  the  new  methods  discussed  in  this  thesis  will 
have  a  strong  indication.  This  chapter  establishes  strengths  and  weaknesses  of  existing 
signal  detection  methods  both  in  the  time  and  frequency  domains.  The  intent  is  not  to 
produce  an  all-inclusive  list  of  current  signal  processing  techniques,  but  rather  to 
compare  new  methods  against  existing  ones.  In  this  way,  we  hope  to  show  that  these  new 
techniques  will  not  replace  current  ones,  but  rather  support  them  as  they  add  to  the 
robustness  of  the  whole  Acoustic  Transient  Detection  System. 

A.  ESTABLISH  A  BENCHMARK  SIGNAL 

In  order  to  meet  the  goal  of  enhancing  transient  signal  detection,  a  comparison 
between  established  methods  and  new  methods  presented  in  this  thesis  must  be 
performed.  To  accomplish  this  evaluation,  a  single  representative  signal  will  be  used  as  a 
benchmark.  The  model  signal  is  the  sound  of  a  fan  blowing  for  six  seconds  sampled  at 
11.025  kHz.  Starting  after  the  third  second,  a  series  of  11  tapping  sounds  (each 
approximately  .1  seconds  in  duration)  will  occur.  These  taps  will  be  part  of  the  analysis 
on  transient  resolution.  The  first  three  seconds  will  be  used  to  study  the  effectiveness  of 
all  transient  signal  detection  techniques  by  varying  the  amplitude,  duration,  and 
frequency  content  of  an  added  transient  signal.  Figure  II- 1  shows  a  block  diagram  of  the 
filtering  process  of  the  signal  to  detect  the  transient  noise.  Figure  II-2  is  a  pictorial 
representation  of  the  benchmark  comparison  setup  based  on  the  filtering  process 
described  in  Figure  II- 1. 
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Figure  II-l  Block  Diagram  of  Whitening  Filter  Process  of  Benchmark  Signal 
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Figure  II-2  Whitened  and  Filtered  Acoustic  Signal  Benchmark  (z[n)) 


1. 


What  Kind  of  Noise? 


It  is  widely  held  that  since  acoustic  ambient  noise  is  comprised  of  many  random 
sources,  it  can  be  assumed  to  be  Gaussian  [2],  This  assumption  is  based  on  the  Central 
Limit  Theorem ,  which  states  that  as  the  number  of  random  variables  in  a  sequence 
becomes  large,  the  probability  distribution  of  the  sum  of  the  sequence  approaches  a 
Gaussian  random  variable  [4],  In  the  case  of  shallow  water,  however,  ambient  ocean 
noise  coming  from  the  Infrasonic  Band  and  the  Low  Sonic  Band  may  dominate,  and 
therefore  may  compromise  the  previous  assumption. 

The  Infrasonic  Band  comprises  frequencies  between  1  Hz  and  20  Hz.  This  is  the 
frequency  band  that  contains  blade  rate  fundamental  frequencies  of  power-driven  vessels. 
The  Low  Sonic  Band  has  a  frequency  range  between  20  Hz  and  200  Hz,  which  is  mostly 
caused  by  distant  shipping,  especially  in  the  frequencies  around  100Hz  [2]. 

At  any  rate,  noise  can  be  defined  simply  as  any  part  a  signal  that  you  don’t  want. 
For  the  specific  case  of  detecting  transient  signals  in  shallow  water,  it  is  a  good 
assumption  that  the  above-mentioned  ambient  noise  bands  dominate  the  noise  model. 

Take  for  example  the  scenario  of  attempting  to  detect  a  small  ship  in  shallow 
water  laying  mines  among  a  fleet  of  fishing  vessels.  It  is  desirable  to  detect  a  certain 
transient  signal  like  a  splash  or  machinery  among  many  acoustic  sources  with  similar 
frequency  content.  Since  we  know  that  mines  are  deployed  in  heavily  traveled 
waterways,  you  can  assume  that  the  background  noise  (defined  as  the  part  of  the  signal 
you  don’t  want)  will  have  a  dominant  low  frequency  content.  The  model  therefore,  will 
not  be  white  noise,  but  colored  noise  since  it  has  uneven  frequency  content. 

The  colored  noise  assumption  is  reflected  in  the  Benchmark  Signal.  Much  like  a 
power-driven  vessel,  a  fan  blowing  produces  a  similar  low  frequency  content  due  to  the 
blade  rate  of  the  fan  as  it  rotates  like  a  propeller. 
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2.  Are  Shallow  Water  Acoustic  Signals  Stationary? 

Like  most  random  processes  occurring  in  nature,  the  background  ocean  ambient 
noise  is  not  stationary.  In  the  example  above,  one  noise  source  out  of  many  may  change 
its  course  or  speed.  So  by  one  slight  change  in  frequency,  a  change  in  the  stochastic 
nature  of  the  ambient  noise  will  take  place.  However,  as  in  many  applications  to  linear 
filtering  techniques,  a  process  can  be  assumed  stationary  over  a  short  period  of  time  [5]. 

The  benchmark  signal  for  this  thesis  assumes  that  the  signal  is  stationary  for  the 
entire  6  seconds.  As  the  results  will  show  (Figure  EE-2),  assuming  the  data  set  is 
stationary  over  this  time  frame  is  sufficient  to  distinguish  the  transient  signal  from  the 
background  noise.  Like  the  Benchmark  Signal,  the  short-time  stationarity  assumption  for 
underwater  acoustics  will  apply. 


B.  STOCHASTIC  SIGNAL  MODEL 

The  first  of  two  current  methods  to  be  discussed  is  an  analysis  in  the  time  domain. 
We  wish  to  transform  the  input  signal  through  a  filter  to  extract  a  transient  from  noise. 
To  do  this,  we  rely  on  the  innovation  process  where  the  output  of  the  filtered  sequence  is 
white  noise. 

It  is  well  known  that  a  process  whose  complex  spectral  density  function  satisfies 
the  Paley- Wiener  Condition  can  be  modeled  as  an  Autoregressive-Moving  Average 
(ARMA)  process  [5],  One  of  the  properties  of  the  signal  model  is  that  it  can  be  reversed 
such  that  it  becomes  a  whitening  filter.  Furthermore,  since  its  transfer  function  Hca  (z)  is 

minimum-phase,  its  inverse  Hc~\z)  will  be  causal  as  well.  Figure  E-3  shows  a  diagram 
of  the  innovation  representation. 
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Figure  II-3  Innovation  Representation  of  a  Random  Process,  (a)  Signal  Model 

(b)  Inverse  Filter  [5] 


1.  Autoregressive  -  Moving  Average  (ARMA)  Model 

Any  regular  stationary  random  process  can  be  represented  as  the  output  of  a  linear 
shift-invariant  filter  driven  by  white  noise  [5].  Starting  with  the  assumption  that  our 
input  signal  x[n]  in  Figure  II-3(b)  is  colored  noise,  we  want  to  design  a  filter  H(z)  such 
that  the  result  will  be  white  noise.  The  innovation  process  will  whiten  the  time  domain 
signal  except  where  the  transients  are  located  because  it  is  not  part  of  the  innovation 
process.  A  second  look  at  Figures  II- 1  and  II-2  will  show  this  process. 

We  use  an  autoregressive  moving  average  (ARMA)  model  in  order  to  represent  a 
random  process.  The  ARMA  model  involves  both  a  nontrivial  numerator  and 
denominator  and  when  compared  to  the  AR  or  MA  models  alone,  they  often  require 
considerably  fewer  parameters  to  match  a  desired  power  spectral  density  function  [5]. 
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Without  a  priori  knowledge  of  the  autocorrelation  function,  we  estimate  it 
directly  from  the  data.  The  Autoregressive  model  estimation  can  then  be  written  as  a 
weighted  sum  of  prior  values  of  the  observations.  Written  as  an  equation,  we  get 


4«]  =  -axx{n  - 1]  -  a2x[n  -2  a  px[n-P],  (2.1) 


where  the  error  is 

P 

£[n]  =  x[n]  -  4«]  =  ^akx[n-k].  (2.2) 

*=0 


In  equation  2.1,  x[n]  is  the  estimate  of  the  n,h  element.  As  in  many  applications, 
we  want  to  minimize  the  error  by  applying  the  Orthogonality  Principle  which  states  that  a 
vector  of  weighting  coefficients  a  can  be  chosen  to  minimize  the  mean-square  error 

E{|x-*|2}if  the  error  is  orthogonal  to  the  observations.  The  final  result  is  the  least 
squares  form  of  the  Yule-Walker  equations, 
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and  S  is  the  sum  of  the  squared  errors.  It  follows  then  that  an  estimate  for  the 
prediction  error  variance  is 


(2.4) 


where  Ns  is  the  size  of  the  data  set. 

The  observation  matrix  X  is  set  up  such  that  it  assumes  the  observations 
outside  of  Ns  data  points  are  zero,  which  is  analogous  to  performing  a  rectangular 

window  operation  on  the  data  set.  Although  this  might  not  be  the  best  method  to  use  in 
terms  of  pole  placement  of  the  vector  a,  the  matrix  structure  has  a  nice  property  since 
(X*rX)  is  Toeplitz.  This  property  means  that  the  matrix  will  be  positive  definite,  which 
guarantees  that  the  whitening  filter  is  minimum  phase  and  therefore  is  a  causal  stable 
system  with  a  causal  stable  inverse.  This  property  is  important  because  we  are  interested 
in  the  innovation  process.  Finally,  using  the  autocorrelation  method  of  determining  the 
Toeplitz  correlation  matrix  (X*rX) ,  fast  algorithms  for  determining  the  elements  of  a 
can  be  used  such  as  the  Levinson  Recursion  [5].  From  the  AR  coefficients  found  from 
equation  2.4,  the  Moving  Average  coefficients  b  can  then  be  computed  by  any  number  of 
algorithms  such  as  the  Basic  Prony  Method  [5].  The  ARMA  routine  used  in  this  thesis 
uses  an  algorithm  to  determine  an  optimum  number  of  poles  and  zeros,  P  and  Q, 
respectively. 

As  an  example,  take  the  Benchmark  Signal  and  find  the  AR  filter  coefficients  a 
and  the  MA  coefficients  b  where 


■  i  ■ 

'  1  ' 

a  = 

a, 

and  b  = 

“p. 

h. 

ll 


The  causal  ARMA  filter  is  then 


Harma(z~1)  = 


b0+bxz~  +b2z  —  +  bQz 


1  +  a{z~l  +  a2z~2  +---  +  apz 


-P  ’ 


(2.5) 


To  realize  the  time  domain  solution  to  the  problem  described  in  Figure  II-3(a),  we 


get 


x[n]  =  ^bqw[n-q]-^ apx[n  - p] . 

<7=0  p- 1 


(2.6) 


Inverting  the  model  as  in  Figure  H-3(b)  the  causal  inverse  ARMA  filter  is  simply 


,-p 


h  1+a,z~  +’"+a^~ 

ARMA  b0+b1z-1+b2z~2+-  +  bQz-Q' 


(2.7) 


The  time  domain  realization  to  the  innovation  process  is  finally, 


w[h]  =  Zj apx[n - p ] ~^bqw[n - q] , 

p=l  9= 0 


(2.8) 


Running  the  Benchmark  Signal  through  the  inverse  ARMA  filter  (as  depicted  in 
the  model  in  Figure  II- 1)  results  in  the  innovation  process,  which  is  white  except  where 
the  transients  occur  in  time.  This  is  because  the  ARMA  model  seeks  only  to  model  the 
colored  noise  as  the  innovation  process.  Since  the  colored  noise  input  has  a  stationary 
localized  frequency  content,  the  ARMA  model  can  do  a  good  job  at  placing  poles  to 
model  the  signal.  All  other  signal  content  that  is  statistically  dissimilar  to  the  colored 
noise  is  rejected  and  appears  as  a  “spike”. 

As  a  check,  in  listening  to  the  fan  blow,  a  low  frequency  “hum”  is  distinguishable 
and  the  1 1  tapping  sounds  are  also  recognizable.  After  the  innovation  process,  the 


12 


background  noise  sounds  like  a  heavy  rain  hitting  the  ground.  When  plotted  the 
transients,  indistinguishable  in  the  time  domain  before  the  innovation  process,  are  now 
quite  visible  against  the  white  noise.  Figure  II-4  shows  the  Benchmark  Signal  before  and 
after  the  innovation  process. 
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Figure  II-4  Innovation  Process  of  Benchmark  Signal  (a)  Input  (b)  Output 
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2.  Strengths  of  the  ARMA  Model 

Since  the  innovation  process  is  based  on  a  statistical  model,  the  frequency  content 
of  the  signal  does  not  affect  its  detection,  even  if  the  signal  contains  the  same  frequencies 
as  the  noise.  The  ARMA  model  provides  a  representation  of  the  noise  and  it  can  detect 
signals  that  are  not  compatible  with  the  model.  The  ARMA  method  picked  up  on  all  the 
transients  and  the  height  of  the  spikes  of  the  output  signal  is  proportional  to  the 
“loudness”  of  the  audible  signal,  which  might  be  an  attractive  feature  to  compare  the 
magnitude  of  two  signals.  This  technique  seems  to  perform  well  on  very  short  transients 
and  it  can  resolve  signals  close  to  each  other.  Notice,  for  example,  the  cluster  of  signals 
just  before  the  third  second  in  Figure  II-4  (b),  where  the  ARMA  model  was  actually  able 
to  distinguish  between  three  transients  less  than  one-tenth  of  a  second  apart  from  each 
other. 


3.  Weaknesses  of  the  ARMA  Model 

The  reverse  ARMA  (innovation)  process  is  a  procedure  that  produces  white  noise 
from  a  colored  noise  input.  This  process,  as  we  shall  see,  is  very  effective  in  detecting 
transient  signals  that  are  not  part  of  the  white  noise  model,  appearing  as  spikes  above  the 
noise  threshold  as  seen  in  Figure  II-4  (b).  However  this  method,  although  performing 
well  on  both  narrowband  and  broadband  signals,  has  its  limitations  when  the  power  of  the 
transient  signal  is  at  or  below  the  modeled  noise  power.  As  we  shall  see,  the  fundamental 
idea  of  the  methods  presented  in  this  thesis  is  to  lower  the  white  noise  power  of  the 
innovation  process  in  order  to  extract  transients  with  smaller  signal  power. 
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C.  TIME-FREQUENCY  ANALYSIS 


The  second  method  to  be  discussed  introduces  a  technique  that  combines  both 
time  and  frequency  analysis.  The  process  is  like  a  sliding  window  in  time  taking  a 
snapshot  of  the  frequency  content  of  a  signal,  assuming  that  the  statistical  properties  do 
not  change  within  the  time  duration.  To  view  the  results,  frequency  is  plotted  as  a 
function  of  time  and  is  called  the  spectrogram.  For  USW,  a  similar  technique  is  used 
called  the  “lofargram,”  which  is  currently  the  standard  method  for  sonar  signal 
processing  [5]. 


1.  Spectrogram 

The  Spectrogram  is  based  on  the  discrete-time  Short-Time  Fourier  Transform, 
which  is  an  extension  to  the  basic  Discrete-Time  Fourier  Transform  and  is  defined  as 


X{n,(Q)=  ^jx(m)w(n-m)e  J0m  ,  (2.9) 

m=-~ 


where 

w(n  -  m)  is  the  sliding  window  function  and 

OO 

X(n,0))=  ^ x(m)e~iean  is  the  DTFT. 

m=-°° 


Completely  analogous  to  the  DTFT,  the  discrete-time  STFT  is  continuous  in 
frequency  so  for  digital  processing,  a  method  to  discretize  frequency  content  like  the  FFT 
can  be  implemented  to  make  this  method  computationally  practical.  Figure  II-5  shows  a 
diagram  of  the  STFT  process.  Along  the  horizontal  axis  is  a  sliding  window  process  as  a 
function  of  time.  At  each  Tlast  +  At ,  a  “snapshot”  of  the  frequency  content  of  the  signal 

at  that  instant  in  time  is  taken.  Another  way  to  visualize  the  process  is  to  look  at  the 
frequency  domain.  Along  the  vertical  axis  is  a  set  of  filter  banks,  each  of  equal 
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bandwidth  A/ .  The  signal  is  passed  through  the  filter  bank  at  each  Tlast  +  Ar  and  the 

spectral  content  within  the  particular  bandpass  filter  is  calculated.  The  latter  pictorial 
representation  of  the  STFT  is  particularly  useful  in  drawing  an  analogy  to  the  Continuous 
Wavelet  Transform  discussed  next. 


Figure  II-5  Time-Frequency  Plane  Corresponding  to  the  Short-Time  Fourier 

Transform  [7] 
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The  (discrete)  STFT  on  the  Benchmark  signal  is  provided  in  Figure  H-6.  The 
darker  colors  on  the  spectrogram  represent  frequencies  with  high  power.  Notice  the  low 
frequency  content  below  300  Hz  is  dominant  throughout  the  entire  6-second  sample. 
Constant  frequencies  at  600  Hz  and  900  Hz  are  also  noticeable.  The  most  revealing 
information  on  the  spectrogram  is  the  vertical  lines  in  the  latter  half  of  the  signal 
indicating  the  presence  of  transients.  What’s  more,  the  discrete  STFT  can  provide 
information  on  the  type  of  transient  it  is.  In  this  case,  the  transients  are  broadband 
signals.  In  addition  to  the  standard  Benchmark  Signal,  three  narrowband  transients  were 
added  to  the  first  second  to  study  the  effects  of  narrowband  transient  resolution.  The 
spectrogram  located  two  of  the  three  signals,  one  at  around  5  kHz  ±  500  Hz  and  one 
around  2  kHz  ±  500  kHz.  The  third  narrowband  transient,  however,  was  not  detected 
because  it  was  masked  by  the  low  frequency  noise. 


17 


Figure  11*6  Spectrogram  of  Benchmark  Signal 


2.  Strengths  of  the  Spectrogram 

Among  the  many  strengths  of  the  spectrogram  is  its  colorful  visual  representation. 
On  a  single  plot,  it  provides  information  in  both  the  time  and  frequency  domains.  This 
representation  therefore  gives  not  only  the  time  at  which  the  transient  signal  occurred,  but 
also  its  frequency  content,  which  may  be  vital  in  classifying  the  source.  For  applications 
in  sonar  signal  processing,  the  spectrogram  is  commonly  referred  to  as  the  “lofargram”, 
which  provides  a  useful  “fingerprint”  mapping  of  unique  acoustic  sources  [5].  From  this 
display  alone,  useful  information  such  as  detecting  dominant  narrowband  tonals, 
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harmonics  and  blade  rates  can  be  seen.  Seeing  the  frequency  content  of  a  contact  change 
over  time  can  also  indicate  a  dynamic  event  such  as  a  change  in  course  or  speed. 

The  spectrogram  works  best  on  broadband  transient  signals  because  there  is  a 
larger  contrast  of  darker  colors  in  a  vertical  line  than  there  is  in  the  single  point 
representation  of  a  narrowband  signal.  In  addition,  we  have  shown  that  it  is  possible  to 
have  a  narrowband  transient  signal  masked  by  the  frequency  content  of  the  background 
noise.  Unless  the  background  noise  is  completely  white,  at  least  some  of  the  frequency 
content  of  a  broadband  noise  is  sure  to  appear. 

Compared  to  the  ARMA  model,  the  STFT  only  assumes  stationarity  over  the 
duration  of  the  sliding  window  only.  Recall  that  the  ARMA  process  attempts  to  model 
the  signal  through  the  statistics  of  the  entire  6-second  data  set.  If  it  were  to  assume  a 
shorter  data  set  to  compute  the  autocorrelation  matrix,  the  whitening  model  will  be 
different  and  may  or  may  not  provide  better  results.  A  further  discussion  on  the  window 
size  will  be  discussed  in  the  final  chapter. 

3.  Weaknesses  of  the  Spectrogram 

Like  the  innovation  process,  the  discrete  STFT  must  assume  stationarity  for  the 
duration  of  the  windowing  function.  There  exist  natural  and  man-rmade  signals  for  which 
its  spectral  content  is  changing  so  rapidly  that  finding  a  proper  windowing  function  is 
difficult  because  there  may  not  be  any  obtainable  time  interval  for  which  the  signal  is 
stationary  [6],  The  result  is  a  transient  signal  may  go  completely  unnoticed. 

Even  if  it  is  a  good  assumption  that  a  short  time  duration  is  stationary,  there  is  a 
tradeoff  to  the  resolution  in  the  frequency  domain  and  vice  versa.  This  is  called  the 
Uncertainty  Principle  or  Heisenberg’s  Inequality,  which  states  that  resolution  in  time  and 
frequency  is  bound  by  the  inequality 


At-Af 


(2.10) 


with  At  and  Af  the  duration  of  the  signal  in  time  and  frequency  respectively. 


19 


If  localization  in  time  is  necessary,  At  can  be  made  small  but  according  to  the 
inequality  (2.10)  A/ will  become  greater  and  therefore  the  spectrogram  will  lose 
frequency  resolution  proportionately. 

To  further  illustrate  this  principle,  the  two  visible  narrowband  transient  signals  in 
Figure  TE-5  are  about  one-tenth  of  a  second  in  duration  and  contain  only  one  frequency 
each,  5  kHz  and  2  kHz  respectively.  From  the  figure,  both  transients  are  stretched 
vertically,  blurring  the  actual  frequency  spectrum  while  the  time  duration  at  which  they 
occur  remains  short.  If  a  broad  window  were  used  in  the  discrete  STFT,  the  frequency 
would  resolve  to  a  single  point  while  the  time  of  the  signal  would  be  stretched  and 
distorted. 

The  spectrogram  merely  shows  the  frequency  of  a  signal  during  a  time  duration. 
Since  the  spectrogram  is  not  based  on  statistical  modeling,  it  cannot  “extract”  a  signal 
from  the  noise  and  therefore  is  unable  to  distinguish  between  signal  and  noise  if  they  both 
have  the  same  frequency  content.  This  is  particularly  important  in  cases  where  an 
acoustic  transient  signal  is  masked  by  the  colored  noise.  Consider  the  Benchmark  Signal 
spectrogram  in  Figure  II-5.  Also  included  in  the  signal  is  a  sine  wave  at  a  frequency  of 
100  Hz.  Since  the  power  of  the  colored  noise  is  highly  concentrated  at  frequencies  below 
300  Hz,  the  transient  signal  goes  undetected.  Unmasking  transient  signals  in  this  region 
become  important  for  sonar  operations  because  much  underwater  activity  may  take  place. 

Finally,  while  the  innovation  process  forces  a  signal  into  a  white  noise  process  in 
order  to  detect  transients,  the  spectrogram  takes  advantage  of  the  signal’s  frequency 
content  over  time.  If  the  input  signal  is  white  noise,  nothing  is  distinguishable  since  the 
frequency  content  is  spread  evenly  over  all  frequencies  for  all  time. 
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III.  WAVELET  ANALYSIS 


Since  we  know  that  the  STFT  is  resolution-limited  in  both  time  and  frequency 
simultaneously,  we  wish  to  investigate  other  means  to  localize  transient  signals  in  time 
while  somehow  still  being  bounded  by  the  uncertainty  principle.  In  this  section,  we  start 
with  a  discussion  of  the  Continuous  Wavelet  Transform  (CWT)  as  it  relates  to  the  Short- 
Time  Fourier  Transform.  From  there,  we  discuss  the  development  of  the  Discrete 
Wavelet  Transform,  its  usefulness  with  filter  banks,  and  the  idea  of  multiresolution 
analysis,  which  is  the  basis  of  the  methods  of  detecting  transient  signals. 


A.  CONTINUOUS  WAVELET  TRANSFORM 

To  describe  the  Continuous  Wavelet  Transform,  we  start  with  the  equation  for  the 
Short-Time  Fourier  Transform 

+oo 

X(f)  =  J x(t)w\t-z)e-2jmdt .  (3.1) 

^O© 

The  integral  can  be  viewed  as  an  inner  product  operating  on  data  inside  the 
windowing  function,  which  measures  the  “similarity”  between  the  function  x(?)  and  a 
sinusoidal  basis  function.  The  result  turns  out  to  be  the  energy  level  at  each  specific 
frequency.  STFT  analysis  is  useful  for  narrowband  periodic  signals  because  of  its  ability 
to  detect  strong  similarity  between  the  input  signal  and  the  sinusoidal  basis  function  at  a 
particular  frequency,  but  it  has  its  limitations  to  the  larger  family  of  non-stationary, 
broadband  signals.  Based  on  Fourier  Transform  restrictions,  it  is  desirable  to  develop  a 
different  transform  that  can  overcome  the  resolution  limitation. 

In  order  to  investigate  a  different  set  of  transforms  that  exhibit  these  properties, 
we  need  to  look  at  other  possible  basis  functions.  To  begin,  assume  that  the  width  of  the 
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frequency  resolution  of  the  transform  A/  is  no  longer  constant  but  rather  the  ratio  of  A/  to 
the  frequency  / ,  resulting  in 


(3-2) 


This  constant  ratio  means  that  as  the  analysis  of  frequency  content  of  a  signal 
increases,  so  does  the  width  of  the  frequency  resolution.  Looking  back  at  the  diagram  of 
the  time-frequency  plane  of  the  STFT  (Figure  D-5),  the  vertical  axis  represents  a  set  of 
constant-width  bandpass  filters  in  a  filter  bank.  In  the  case  of  the  CWT,  the  constant- 
width  frequency  resolution  is  replaced  by  a  filter  bank  of  constant  relative  bandwidth 
(called  “constant-Q”  analysis).  This  process,  in  essence,  allows  the  resolution  in  At  and 
A/  to  vary  in  the  time-frequency  plane  in  order  to  obtain  a  multiresolution  analysis, 
which  will  be  discussed  in  greater  detail  in  the  next  section.  Another  way  to  look  at  the 
filtering  process  is  that  the  filter  bank  used  is  not  linearly  spaced  along  the  frequency  axis 
(as  for  the  STFT  case)  but  spaced  evenly  over  a  logarithmic  scale  [7].  Figure  III.  1  shows 
this  process.  Now  substituting  (3.2)  into  the  Heisenburg  Inequality  (2.10)  gives 


At  = 


1 

4 7tfc  ‘ 


(3.3) 


By  looking  at  equations  (3.2)  and  (3.3),  one  can  see  that  as  the  frequency  becomes 
greater,  the  resolution  in  time  decreases  and  the  resolution  in  frequency  increases. 
Conversely,  as  the  frequency  becomes  smaller,  the  resolution  in  time  becomes  larger  and 
the  resolution  in  frequency  becomes  smaller.  This  is  a  fundamental  property  of  the  CWT. 
The  time  resolution  of  the  CWT  becomes  arbitrarily  good  at  high  frequencies  while  the 
frequency  resolution  becomes  arbitrarily  good  at  low  frequencies. 
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Figure  III-l  Filter  Division  in  the  Frequency  Domain,  (a)  Uniform  Coverage  (STFT) 

(b)  Logarithmic  Coverage  (CWT)  [7] 


The  definition  of  the  Continuous  Wavelet  Transform  is  thus 
CWTX  (t, a)  =  J x(t)h\,r  ( t)dt , 


(3.4) 


where 


't-T' j 

—  ll 

l\ 

,  J 

(3.5) 
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is  the  basis  function  called  a  wavelet. 

A  description  of  a  wavelet  function  will  be  addressed  in  the  next  section.  For 
now,  some  of  the  properties  of  the  transform  will  be  discussed.  From  the  definition,  it 

can  be  seen  that  a  is  the  “scaling”  factor  described  in  equation  (3.2)  and  l/^ja[  is  used 

for  energy  normalization.  The  numerator  of  the  wavelet  function  shows  the  time  -  shift 
property  of  the  sliding  window  function. 

By  way  of  analogy,  the  STFT  is  plotted  as  frequency  vs.  time  as  a  spectrogram. 
The  CWT  is  a  function  of  time  and  scale  and  therefore  will  be  plotted  as  scale  vs.  time  as 
a  “scalogram”.  Figure  113-1  is  a  comparison  of  a  CWT  to  a  STFT  of  the  same  signal. 
Notice  that  the  CWT  of  the  delta  function  converges  to  a  single  point  at  t  =  t0  and  the 

energy  is  spread  out  with  increasing  scale.  The  STFT  of  the  delta  function  shows  the 
resolution  limitation  in  the  time  domain  with  equal  energy  at  all  frequencies.  For  the 
sinusoidal  signal  of  three  frequencies  at  /0 ,  2  f0  and  4  f0 ,  the  STFT  shows  the  constant 

frequency  in  time,  bound  by  the  resolution  limitation  in  the  frequency  domain.  For  the 
CWT,  however,  we  see  the  effect  of  the  logarithmic  filter  bank  where  the  low  frequency 
resolution  is  narrowest  and  continues  to  double  in  size  as  the  frequency  increases.  Figure 
m-2  shows  clearly  the  principle  that  the  time  resolution  of  the  CWT  becomes  arbitrarily 
good  at  high  scale  while  the  frequency  resolution  becomes  arbitrarily  good  at  low  scale. 
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Figure  III-2  Comparison 


CWT  to  the  STFT  [7] 


B.  DISCRETE  WAVELET  TRANSFORM 


At  this  point  a  discussion  of  the  properties  of  the  Continuous  Wavelet  Transform 
has  been  made  as  it  relates  to  the  STFT  and  the  benefits  in  localization  in  both  time  and 
scale  have  been  presented.  The  CWT  also  reduces  the  number  of  coefficients  necessary 
to  define  the  signal.  During  the  discussion,  ideas  of  orthogonality,  energy  distribution, 
logarithmic  filter  banks,  scaling  functions,  and  wavelet  basis  functions  emerged  as 
important  properties  of  the  CWT.  This  section  will  continue  to  build  on  these  concepts  as 
they  apply  to  the  Discrete  Wavelet  Transform  (DWT). 

For  reasons  of  clear  analysis  and  better  understanding,  we  begin  with  the 
definition  of  a  linear  decomposition.  A  real-valued  function  can  be  expressed  as  a  linear 
decomposition  such  that 


/(0  =  (3-6) 

i 

Parameter  a;  is  the  real-valued  expansion  coefficients,  and  is  a  set  of  real¬ 
valued  functions  of  t  called  the  expansion  set.  When  the  expansion  is  unique,  the  set  of 
functions  y/t(t)  is  called  a  basis  and  it  is  used  to  represent  f(t) .  If  the  basis  is 
orthogonal,  i.e. 


{Vk>Vi)  =  \ Vk ( t)Wi 0 t)dt  =  0  ,  k*l  (3.7) 


then  the  coefficients  can  be  determined  by  the  inner  product 

ak={fit),¥k)  =  \f{t)y/k{t)dt.  (3.8) 


For  the  wavelet  expansion,  however,  a  two-parameter  system  is  used  such  that  the 
basis  used  to  represent  f{t)  is  now 
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(3-9) 


*  j 


where 


j  represents  the  integer  scaling  index  of  the  scaling  function, 
k  represents  the  integer  time  translation, 

ilfj  k  ( t )  are  the  wavelet  expansion  functions  that  usually  form  an, 


orthogonal  basis,  and 

aJ  k  are  the  set  of  expansion  coefficients. 


The  wavelet  basis  set  y/j  k  (t)  can  be  defined  in  terms  of  scaled  and  translated 
versions  of  the  Mother  Wavelet  y/(t) ,  i.e. 

y/.k(t)  =  2^(2jt-k).  (3.10) 


Here,  the  term  2 2  is  a  normalizing  factor  as  j  increases.  This  property  described 
in  equation  (3.10)  is  called  the  multiresolution  condition,  which  states  that  if  a  set  of 
signals  can  be  represented  by  a  weighted  sum  of  <p(t  -  k) ,  then  a  larger  set  (including  the 
original)  can  be  represented  by  a  weighted  sum  of  (p{2t-k).  In  other  words,  “if  the 
basic  expansion  signals  are  made  half  as  wide  and  translated  in  steps  half  as  wide,  they 
will  represent  a  larger  class  of  signals  exactly  or  give  a  better  approximation  of  any 
signal”  [9]. 
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1.  Signal  Spaces  and  Multiresolution  Analysis 

To  discuss  the  DWT,  it  is  best  to  start  with  the  idea  of  resolution  to  define  the 
effects  of  changing  scale.  To  that  end,  we  start  with  a  scaling  function  <p(t ) ,  which  then 
will  define  the  wavelet  basis  function  y/(t) .  We  define  a  set  of  scaling  functions  in  terms 
of  integer  translates  of  the  basic  scaling  function  by 

(pk(t)  =  <p(t-k)  keZandtpeL2.  (3.11) 

where  the  L2  space  is  defined  as  the  set  of  all  functions  f(t)  with  finite  energy, 
i.e., 


J|/(0|2<*<°°- 

The  subspace  spanned  by  the  functions  <pk  (?)  is  v0  such  that 

v0  =  Span{<pk  (?),  keZ}.  (3.12) 

That  is  to  say,  the  subspace  v0  is  spanned  by  the  closed  set  of  all  functions  <pk  (?) 
over  all  integer  values  k .  From  equation  (3.12)  we  can  say  that  a  function  fit)  within 
the  subspace  v0  can  be  expressed  as  a  linear  combination  of  the  set  of  functions  <pk  (?) 
such  that 


/(0  =  Xa^(0  for  any  /(?) 6v0.  (3.13) 

k 

From  here,  we  consider  a  two-dimensional  family  of  functions  and  derive  a  similar 
linear  combination  as  in  equation  (3.13).  We  start  with  the  definition  of  a  scaling 
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function  parameterized  by  scaled  and  translated  versions  of  (pit)  similar  to  equation 
(3.10),  resulting  in 

i 

(pjk{t)  =  22(p{2Jt-k) ,  ke  Z  (3.14) 

where 

Vj  =  Span\ff>j  k  {t),  k  e  z}.  (3.15) 

In  the  same  way  v0  is  expressed  as  a  linear  combination  of  the  set  of  functions 
<pk{t)  inequation  (3.13),  fit)evj  can  be  expressed  as 

/(0  =  5X$*2'r  +  *).  (3-16) 

k 

From  here  we  turn  to  a  useful  property  of  the  wavelet  expansion  set  [9].  There  is 
a  basic  constraint  of  multiresolution  analysis,  which  requires  a  “nesting”  of  the  spanned 
spaces  as 

...cv_2  cv.,  cv0  cv,  cv2  C...L2, 

or  simply 

Vj  c  vy+I  for  all  j'eZ.  (3.17) 

By  this  observation,  in  order  to  ensure  that  the  elements  in  v;  are  scaled  versions 
of  the  next  higher  space  vj+1 ,  we  add  the  constraint. 
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mevj 


(3.18) 


<=>  f(2t)evj+l. 

This  is  done  to  ensure  that  analysis  and  synthesis  can  be  performed.  That  is  to 
say,  if  we  conduct  a  wavelet  decomposition  at  a  certain  resolution  (or  scale),  there  exists 
a  reversing  process  whereby  a  signal  can  be  reconstructed  (synthesized)  to  its  original 
state.  Figure  IH-3  shows  a  Venn  diagram  of  the  nested  vector  space  property  of 
multiresolution  analysis. 


This  set  of  nested  vector  spaces  now  meets  the  multiresolution  condition,  which 
states  that  if  a  set  of  signals  can  be  represented  by  a  weighted  sum  of  <Pj  k  (t) ,  then  a 

larger  set  can  be  represented  by  a  weighted  sum  of  <pj+1Jc  ( t )  .  In  order  to  relax  notation, 
we  drop  k  and  assume  that  (pj  (r)  is  dependent  on  the  time  translation.  The  nesting  of  the 
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spans  as  shown  in  Figure  IE-3  and  by  equations  (3.17)  and  (3.18)  is  obtained  by  requiring 
that  <p0(t)  is  not  only  in  the  space  v0  but  also  in  vx ,  which  is  spanned  by  the  function 

(pl  (t) .  From  equation  (3.16),  we  can  see  that  that  the  function  <p0  (t)  within  the  subspace 
v0  can  be  expressed  as  a  linear  combination  of  the  set  of  functions  <px  (t)  such  that 


Po(0  =  5X^- 

n 


(3.19) 


Substituting  <p0(t)  and  <pt  (t)  using  equation  (3.14)  we  obtain  an  equation  in  terms 
of  scaled  and  translated  versions  of  the  scaling  function  as 

<p(t)  =  y£h(n)y[2<p(2t-n),  ne  Z.  (3.20) 

n 

This  equation  shows  a  telescoping  structure  of  the  scaling  function  and  is  called 
the  Multiresolution  Analysis  (MRA)  equation,  where  h(n)  represents  the  weighted 

coefficients  called  the  scaling  function  coefficients. 

We  now  have  a  recursive  scaling  function  spanning  a  subspace,  which  is  also  a 
subset  of  the  next  higher  scale.  There  is  one  more  feature  of  a  signal  that  does  not  use  a 
scaling  function,  but  rather  uses  a  different  set  of  functions,  called  .  This  set  of 
functions  is  defined  as  a  basis  set  that  spans  the  differences  between  the  spaces  spanned 
by  the  various  scales  of  the  scaling  function  [9].  As  we  shall  see,  these  are  the  set  of 
wavelet  functions. 

Although  wavelets  and  their  corresponding  scaling  functions  are  not  required  to 
be  orthogonal,  for  the  purposes  of  this  thesis,  we  will  limit  our  discussion  to  the  set  of 
orthogonal  basis  functions.  The  advantages  of  having  an  orthogonal  basis  are  that  the 
wavelet  expansion  coefficients  are  easy  to  calculate  and  are  governed  by  Parseval’s 
theorem  which  allows  the  partitioning  of  signal  energy  in  the  wavelet  transform  domain, 
like  that  of  the  classic  Fourier  transform  domain.  In  vector  function  space  notation,  we 
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define  as  the  orthogonal  complement  of  vy.  in  vj+1 .  That  is  to  say,  {w;}  is  the  set  of 
all  functions  that  span  vJ+l  and  are  not  included  in  vy .  v}  is  also  be  considered 
uncorrelated  with  Wj  and  the  union  of  the  two  vector  spaces  span  v  .+1 .  Mathematically, 
all  functions  spanned  by  v j  are  orthogonal  to  vv; .  Therefore  the  basis  functions  q>]  k  and 
y/jj  are  orthogonal  as 


{<Pj,k  (0.  ¥j.i  (0)  =  J  (phk  (t)Vjj  (1 t)dt  =  0 . 


(3.21) 


Starting  at  an  arbitrary  scale  spanned  by  the  scaling  functions  in  vy. ,  say  j  =  0, 
we  can  write 


v0  cv,  cv2  c-cl2. 

Now  from  the  orthogonal  relationship  between  v0  and  w0 ,  we  get 

V!  =  V0@w0,  (3.22) 

where  ©  indicates  a  union  operator. 

Finally,  due  to  the  “telescoping”  nature  of  the  subspaces  v;. ,  j  e  Z  and  the 
orthogonality  of  the  corresponding  wy. ,  je  Z ,  all  of  L 2  can  be  expressed  as 

L2  =  v0  ©  w0  ©  Wj  ©  w2  ©  •  •  • .  (3.23) 

Figure  HI-4  shows  the  Venn  diagram  of  the  entire  vector  space  spanned  by  L2. 
Assuming  that  our  arbitrary  starting  point  in  the  infinitely  telescoping  set  of  nested  vector 
spaces  spanned  by  the  scaling  functions  is  v} ,  we  can  see  that  wj  is  orthogonal  to  v. 
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and  a  subset  of  vj+l  as  shown  in  Figure  HI-5.  These  wavelets  reside  in  the  space  spanned 

by  the  next  narrower  scaling  function  and  therefore  the  multiresolution  principle  applies. 
Since  w0  c  v, ,  the  set  of  wavelets  in  w0  can  be  represented  as  a  weighted  sum  of  the 
next  higher  scale  (p{2t)  resulting  in 

y/{t)  =  ^  hy  {n)^2<p{2t  -  n) ,  ne  Z ,  (3.24) 

n 

where  l\(n)  represents  the  weighted  coefficients  that  define  the  wavelet. 

It  can  be  seen  that  in  equation  (3.24)  the  wavelet  function  is  expressed  in  terms  of 
a  linear  combination  of  scaled  and  translated  versions  of  the  scaling  function,  which 
implies  a  relationship  between  the  scaling  function  and  the  wavelet  function.  Indeed,  by 
taking  equations  (3.20)  and  (3.24),  it  can  be  shown  that  h(n )  and  h^in)  are  related  by  [9] 

hl(n)  =  (-l)nh(-n).  (3.25) 
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Figure  III-4  Scaling  Function  and  Wavelet  Vector  Spaces  [8] 
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2.  Implementation  of  the  Discrete  Wavelet  Transform 

We  now  have  defined  a  wavelet  system  in  terms  of  its  scaling  function  and 
Mother  Wavelet.  As  it  turns  out,  it  is  not  necessary  to  actually  know  <p(t)  and  y/(t) .  All 
we  need  are  the  filter  coefficients  to  the  MRA  equation  in  order  to  perform  a  discrete 
wavelet  decomposition.  This  section  introduces  the  Discrete  Wavelet  Transform  (DWT) 
and  how  it  can  be  implemented. 

In  the  last  section  we  described  the  whole  vector  function  space  as  a  union  of  the 
coarsest  scale  with  all  wavelet  subspaces  such  that 


L2  =  v0  ©  w0  ©  w,  ©  w2  ©  •  •  • . 


Based  on  this  fact,  we  can  now  express  a  function  g{t)  as  a  series  expansion  in 
terms  of  the  scaling  function  <pk  (t)  and  the  wavelet  functions  y/jk  (t)  giving  the  general 
form  of  the  equation 

*(')  =  «>  +  £  IX  (*#M  («> ■  (3.26) 

k  k  j=jO 

Here,  jO  is  the  arbitrary  position  for  the  coarsest  scale  whose  space  is  spanned  by 
the  scaling  function  <pj0  k  (t) .  Notice  that  the  inner  summation  on  the  second  term  begins 
at  JO ,  indicating  that  the  wavelet  coefficients  begin  at  the  same  level  as  the  current  scale 
level.  If  the  conditions  of  orthogonality  are  met,  these  wavelet  coefficients  can 
completely  describe  the  function  g(t)  and  are  useful  for  analysis,  denoising,  compression 
and  perfect  reconstruction  of  the  original  signal  [10].  Assuming  that  the  wavelet  system 
is  orthogonal,  the  coefficients  Cj(k)  and  d}  (k)  can  be  found  by  the  inner  product  of  (pj  k 

and  y/jk  respectively. 

Performing  the  inner  product  to  find  Cj  ( k )  we  have 
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Cj  (k)  =  (g(t),<pjtk  (i t )}  =  J  g(t)<pjk  (t)dt 
=  J g(i)2}(p(2jt-k)dt .  (3.27) 

Restating  the  MRA  equation, 

<p(t)  =  h{n)4l(p{2t  -  n ) , 

n 


we  substitute  in  the  next  higher  scale  and  translate  the  function  by  k,  resulting  in 
<p(2Jt-  k )  =  J  Kn)y[2<p(2(2i  t  -  2k  -  n)) . 

n 


By  making  a  variable  change  m  =  2k  +  n ,  we  get 

<p(2jt-k)  =  '£h(m-2k)^<p(2j+1t-k).  (3.28) 

m 


Now  substituting  (3.28)  into  (3.27)  it  can  be  shown  that 


cj(k)  =  ^jh(m-2k)jg(t)2  2  <p(2j+lt-k)dt.  (3.29) 


Finally,  the  integral  expression  is  simply  the  inner  product 

(g(tl<Pj+i,m(t))  =  cj+l(m), 

and  therefore 
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(3.30) 


cj  (k)  =  ^h(m-  2  k)cj+l  (m) . 

m 

A  similar  derivation  for  finding  the  coefficients  d  ■  ( k )  yields 

dj(k)  =  ^hl(m- 2k)cj+l  (m) .  (3.3 1) 

m 

Equations  (3.30)  and  (3.31)  show  that  the  scaling  and  wavelet  coefficients  at  the 
next  lower  scale  j  can  be  computed  by  convolving  the  scaling  coefficients  at  the  current 

scale  y  +  1  with  the  time-reversed  weighting  coefficients  h(-n)  and  1\  (-n) ,  then 
downsampling  the  result.  Figure  III-5  shows  a  diagram  of  this  process. 


Figure  III-6  Two-Stage,  Two-Band  Analysis  Tree  [8] 
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3.  Analysis  in  the  Frequency  Domain 

In  the  following  section,  it  will  be  shown  that  the  coefficients  of  the  Haar  Wavelet 
filter  coefficients  h0  (n)  are  a  lowpass  filtering  operation.  Likewise  the  coefficients  \  (n) 

are  a  highpass  operation  [9].  Together,  they  form  a  Quadrature  Mirror  Filter  Bank 
(QMF)  system.  In  addition,  a  reverse  filtering  operation  can  be  performed  allowing  for 
perfect  reconstruction  of  the  original  signal,  also  known  as  a  synthesis  procedure.  This 
makes  sense  because  for  any  transform,  a  reverse  operation  must  be  obtainable  for  the 
transform  to  be  widely  used. 

In  the  frequency  domain,  the  signal  coming  in  is  processed  through  a  filter  bank 
that  splits  the  frequency  content  into  two  parts.  The  high  frequency  band  of  the  signal  is 
represented  by  the  wavelet  coefficients  also  known  as  the  details.  The  low  frequency 
portion  of  the  signal  is  represented  by  the  scaling  coefficients  also  known  as  the 
approximations.  The  lowpass  portion  can  then  be  processed  through  the  same  filter  bank 
again.  Theoretically,  lowpass  filtering  the  scaling  coefficients  can  be  repeated 
indefinitely,  but  in  practice  is  limited  to  the  number  of  data  points  in  the  signal.  This 
limitation  exists  because  every  time  the  signal  is  passed  through  the  filter  bank  it  gets 
decimated  by  a  factor  of  2,  which  effectively  means  the  number  of  data  points  in  the 
signal  is  cut  in  half. 

The  first  stage  splits  the  frequency  band  in  half,  the  second  stage  splits  the  low 
frequency  portion  into  quarters,  and  so  on.  The  resulting  “analysis  tree”  makes  up  a 
logarithmic  filter  bank  as  seen  in  the  CWT.  A  diagram  of  this  process  can  be  seen  in 
Figure  DI-7  as  it  compares  to  a  constant  width  filter  bank. 
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a) 


Constant  Bandwidth  (STFT  case) 


Figure  III-7  (a)  Constant  Width  Filter  Bank  (b)  Frequency  Bands  for  the  Analysis 


Tree  [8] 


39 


C.  HAAR  WAVELET 

Up  to  this  point  we  have  studied  the  Wavelet  Transform  as  a  linear 
decomposition.  We  now  move  from  the  abstract  to  the  application.  Equations  3.30,  3.31 
and  Figure  HI-6  all  describe  the  process  of  decomposing  a  signal  into  successive  details 
and  approximations.  Now  looking  at  the  decomposition  as  an  application  to  signal 
processing,  we  see  that  it  is  a  series  of  simple  filtering  and  downsampling  operations. 
For  our  particular  application,  the  actual  numerical  values  of  the  decomposition 
coefficients  are  not  important,  only  their  relative  values  and  so  the  coefficients  can  be 
normalized.  By  viewing  the  DWT  as  a  simple  signal  processing  procedure,  we  have 
essentially  found  a  way  to  split  and  downsample  a  signal  into  mutually  independent 
signals.  As  we  shall  see,  this  procedure  provides  useful  properties  in  detecting  transient 
signals  by  removing  background  noise. 


1.  Haar  Wavelet  Parameters 

We  are  now  ready  to  discuss  the  actual  wavelet  filter  coefficients  h0(n )  and 
(n) .  There  are  many  wavelet  basis  sets  to  choose  from  and  a  discussion  on  choosing 
specific  types  will  be  discussed  in  the  final  chapter.  For  the  purposes  of  this  thesis,  we 
will  choose  the  Haar  Wavelet  system  and  derive  all  transient  detection  methods  and 
conclusions  based  on  this  set.  Admittedly,  this  is  the  simplest  of  all  wavelet  basis  sets  but 
the  results  will  show  a  marked  improvement  over  current  detection  schemes.  In  addition, 
it  is  more  intuitive  to  describe  the  filtering  processes  in  the  upcoming  chapters  using  Haar 
wavelet  coefficients.  Once  the  transient  detection  methods  have  been  discussed,  a  more 
general  case  can  be  easily  substituted. 

The  Haar  Wavelet  coefficients  are  defined  as 


/*(«)  =  {/z(0),fc(l)}= 


V2  V2 


(3.32) 
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for  the  lowpass  filter  and 


A, 00  =  {*,(0)A<1)}=U.  -j=j  (3.33) 

I 

for  the  highpass  filter. 

A  more  detailed  discussion  explaining  how  the  Haar  Wavelet  meets  the  necessary 
conditions  of  a  Wavelet  system  is  presented  in  Appendix  A.  To  find  the  shape  of  the 
scaling  function,  recall  the  MRA  equation  (3.20) 

(p(t)  =  ^h(n)42(p(2t-n). 

n 

Substituting  the  scaling  function  coefficients,  the  scaling  function  becomes 

(p{t)  =  (p{2t)  +  <p(2t  -  n ) .  (3.34) 

Equation  (3.34)  clearly  shows  that  the  Haar  Scaling  function  is  a  unit-width,  unit- 
height  pulse  function  and  is  equal  to  the  sum  of  two  scaled  and  shifted  versions  of  itself. 

To  find  the  shape  of  the  Haar  Wavelet,  we  use  the  wavelet  equation  (3.24) 

W(t)  =  ^hl(n)^<p(2t-n). 

n 

Substituting  the  Haar  Wavelet  function  coefficients,  the  above  equation  becomes 

y/{t)  =  <p(2t)-<p(2t-\) .  (3.35) 

In  the  same  manner  as  the  scaling  function,  the  shape  of  the  Haar  Wavelet  is  the 
sum  of  two  scaled  and  translated  versions  of  (pit) ,  except  now  the  coefficient  on  the 
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For  the  analysis  of  both  transient  detection  procedures,  it  is  vital  to  understand  the 
effect  of  the  wavelet  decomposition  filtering  process  in  the  frequency  domain.  We  start 
out  by  viewing  the  Haar  Scaling  filter  coefficients  in  the  z  -  domain  as 


(3.36) 


It  is  important  to  note  here  that  the  frequency  response  of  the  scaling  coefficients 
h(n )  (also  known  as  h0  (n) )  is  that  of  a.  lowpass  filter  with  its  cutoff  frequency  at  the 
Nyquist  frequency  rate  (H(n)  =  0).  It  can  also  be  shown  that  the  frequency  response  of 
the  wavelet  coefficients  /ij  (n)  is  that  of  a  highpass  filter  with  its  maximum  at 
H(n)  =  V2  and  its  zero  magnitude  at  H( 0) .  This  is  the  property  that  allows  the  input 
signal  to  be  filtered  into  a  “constant-Q”  filter  bank  resulting  in  equal  highpass  and 
lowpass  portions  preceding  every  downsample  as  described  in  Figure  ni-7.  Figure  III-9 
shows  the  frequency  response  of  the  PR  filter  bank  used  for  the  Haar  Wavelet  system. 


Figure  III-9  Frequency  response  of  Haar  Scaling  and  Wavelet  Functions 


2.  Haar  Wavelet  Summary 

For  the  DWT,  the  Haar  Wavelet  and  Scaling  functions  themselves  are  only 
important  in  terms  of  describing  the  wavelet  decomposition  process  heuristically.  The 
only  parameters  needed  to  apply  wavelet  analysis  are  the  scaling  coefficients  to  the 
Multiresolution  Analysis  equation.  Table  III-l  shows  a  summary  of  the  Haar  Wavelet 
system. 

Indeed,  the  scaling  coefficients  h0  (ri)  and  the  wavelet  coefficients  /i,  (n)  work 
together  to  form  a  perfect  reconstruction  filter  bank.  When  a  signal  is  processed  through 
the  FIR  filter  \  in) ,  it  is  being  passed  through  a  highpass  filter.  When  the  signal  is 
downsampled,  the  set  of  values  is  the  wavelet  decomposition  coefficients  called  the 
“details”.  Likewise  when  the  same  signal  passes  through  the  FIR  filter  h0(n),  it  is  being 

processed  through  a  lowpass  filter.  After  the  signal  is  downsampled,  the  set  of  values  is 
the  scaling  decomposition  coefficients  called  “approximations”.  The  approximations  can 
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then  be  decomposed  again  using  the  same  filter  bank  and  downsampling  indefinitely 
(theoretically)  to  allow  for  infinite  (theoretically)  resolution  in  both  time  and  scale 
(Figure  HI-6). 

Finally,  all  details  at  different  scale  and  the  approximation  at  the  smallest  scale 
are  orthogonal  to  each  other.  This  property,  as  we  shall  see,  implies  that  all  signals  at  the 
output  of  all  filter  bank-downsample  phases  are  uncorrelated  to  each  other  and  is  the 
basis  for  the  first  transient  detection  method. 


NAME  FUNCTION  SHAPE  COEFFICIENTS  FREQUENCY 

RESPONSE 


Table  II- 1  Haar  Wavelet  System  Summary  [15] 
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IV.  METHOD  Is  INNOVATION  -  DWT  PROCESS 


With  the  desirable  properties  of  the  DWT,  we  focus  on  the  first  method  to 
enhance  transient  detection.  This  section  will  present  a  signal  detection  design  based  on 
decomposing  the  innovation  process  through  the  Haar  Wavelet  system.  The  desired 
result  will  be  to  remove  noise  through  the  filtering  process  while  keeping  the  signal 
intact. 


A.  METHOD  1  DERIVATION 

As  mentioned  in  the  previous  section,  the  first  method  takes  advantage  of  the 
orthogonality  of  subspaces  in  wavelet  decomposition.  To  say  that  the  details  and 
approximations  at  each  stage  are  orthogonal  implies  that  they  are  uncorrelated  to  one 
another.  As  a  proof,  consider  the  whitened  input  signal  processed  through  one  stage  of 
the  Haar  Wavelet  system  as  in  Figure  IV- 1.  Note  that  e^il]  and  eHP[l]  are  using  the 
same  index  l  because  they  are  produced  from  the  same  input  and  tire  decimated  by  the 
same  value. 
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Figure  IV-1  Wavelet  Decomposition  on  Whitened  Signal 


From  the  previous  figure,  we  can  express  the  approximations  eLP[l]  and  the 
details  eHP[l ]  in  terms  of  the  input  signal  e[n]  such  that 


eLP  [*] =  ~{e\.^\  +  e[2l  - 1])  and 

(4.1) 

eHp[l]  =  ^(e[2l]-e[2l-l]). 

(4.2) 

To  determine  the  correlation  between  the  two  signals,  we  perform  the  cross 
correlation  operation,  i.e. 

E{.eLP  M  ‘  eHP  =  ~  E{{e\2l]  +  e[2l  —  1])  •  {e[2/]  —  e[2 1  —  1])}  (4.3) 
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Carrying  through  the  multiplication,  the  expression  yields 


E{eu,m-e„M=\AWlf)-\E{Wl-liV°  <4-4) 


Since  the  signal  e[n]  is  wide  sense  stationary,  2?{|e[2Z]|2}  is  equal  to 

e{|^[2/  — 1][ 2 }  and  equation  (4.4)  is  zero  and  hence,  e^ll]  and  eHP[l]  are  uncorrelated. 

The  idea  of  Method  1  is  to  perform  a  DWT  on  the  whitened  signal.  This  setup 
passes  a  white  noise  process  through  the  wavelet  decomposition  where  we  expect  white 
noise  at  each  resolution  stage.  Based  on  the  properties  of  the  Discrete  Wavelet 
Transform  and  the  proof  at  the  beginning  of  the  chapter,  all  the  details  and  the 
approximations  at  the  same  multiresolution  stage  are  uncorrelated.  The  description  here 
is  a  denoising  process  where  we  attempt  to  lower  the  noise  power  by  successive  filtering 
and  downsampling  procedures. 

B.  METHOD  1  SUMMARY 

Method  1  uses  the  Haar  Scaling  and  Wavelet  filter  coefficients  for  the 
multiresolution  filter  bank  system  and  is  shown  in  Figure  IV-2.  Keep  in  mind  this 
method  assumes  that  the  DWT  is  operating  on  a  white  noise  signal,  which  means  that  the 
use  of  the  ARMA  model  is  used  prior  to  the  first  stage  of  the  filter  bank.  If  the 
innovation  process  detects  a  transient  signal,  note  that  there  is  no  need  for  this  algorithm. 
What  this  method  attempts  to  do  is  go  beyond  the  detection  capabilities  of  the  innovation 
process  by  removing  excess  noise  in  order  to  enhance  the  robustness  of  the  detection 
process. 
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12 


>erp[k] 


Figure  IV-2  Method  1:  Innovation  -  Discrete  Wavelet  Transform  Process 
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V. 


METHOD  2:  DWT  -  INNOVATION  PROCESS 


We  now  take  a  different  approach  to  the  same  problem  using  a  derivative  of  the 
Haar  Wavelet  and  its  multiresolution  capability.  As  shown  in  Figure  IV- 1,  Method  1 
took  the  wavelet  decomposition  after  the  innovation  process.  Now  we  will  investigate  a 
method  where  we  perform  the  multiresolution  filter  bank  process  before  the  innovation 
process  as  in  Figure  V-l.  The  idea  of  this  method  is  to  investigate  the  effectiveness  of 
another  denoising  scheme. 

Method  2  takes  a  signal  and  processes  it  through  the  wavelet  decomposition.  By 
the  properties  of  the  DWT,  the  signals  at  each  channel  will  be  mutually  uncorrelated.  At 
this  point,  the  decomposed  signal  is  still  not  whitened.  In  order  to  whiten  the  input,  there 
must  be  an  independent  ARMA  filter  for  each  channel  because  the  same  set  of  ARMA 
filter  coefficients  is  only  useful  for  one  channel.  Fortunately,  it  turns  out  that  just  as  there 
is  a  relationship  between  the  filters  in  the  wavelet  filter  bank,  there  is  a  corresponding 
relationship  between  ARMA  filters  in  the  innovation  filter  bank  and  is  dependent  on  the 
Haar  Wavelet  filter  coefficients. 

This  chapter  will  discuss  a  procedure  to  determine  a  set  of  whitening  filter 
coefficients  from  the  original  set  of  ARMA  filter  coefficients.  The  first  section  starts  out 
with  modified  Haar  Wavelet  filter  bank  and  begins  the  procedure  to  determine  all  the 
ARMA  filter  coefficients. 
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A.  METHOD  2  DERIVATION 


Figure  V-l  shows  a  block  diagram  of  the  multiresolution  innovation  process. 
First,  note  that  the  multiresolution  filter  bank  is  a  generic  highpass/lowpass  filter  pair  that 
are  orthogonal  to  one  another.  Instead  of  the  Haar  Wavelet,  the  filter  coefficients  are 

multiplied  by  a  scalar  value  of  -7=  so  that 


h0(n)  = 


s 


lo,Haar 


(n)  = 


(5-1) 


and 


hl(n)  = 


\,Haar  (n)  ~ 


2 


(5.2) 


The  new  set  of  highpass  and  lowpass  filters  can  be  viewed  as  an  averaging 
(lowpass)  filter  and  a  difference  averaging  (highpass)  filter.  The  innovation  at  the  output 
of  the  whitening  filter  bank  are  orthogonal  to  each  other  where  the  “tildes”  represent  the 
highpass  filter  channels,  or  in  terms  of  the  DWT,  the  details.  The  subscripts  represent  the 
relative  sampling  frequency  in  the  sense  that  ek  is  decimated  by  a  factor  of  2* .  The 

dashed  lines  in  Figure  V-l  represent  a  set  of  cascaded  lowpass  filters  that  are 
intermediate  steps  for  further  decomposition,  yet  can  also  provide  useful  information. 
Although  the  decomposition  can  go  on  indefinitely,  there  comes  a  point  where  a  short 
transient  signal  is  averaged  excessively  with  the  background  noise.  By  experimentation, 
three  levels  of  resolution  were  found  to  be  the  most  useful  in  this  thesis. 
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It  is  also  important  to  note  here  that  the  wavelet  and  ARMA  filter  coefficients  are 
interrelated,  and  the  model  is  designed  exclusively  for  the  Haar  wavelet  and  its 
derivative.  A  more  general  form  of  Method  2  is  possible,  allowing  the  use  of  any  wavelet 
system,  but  is  a  subject  for  further  research. 

Consider  the  ARMA  process  without  multiresolution  as  in  chapter  2.  For  notation 
simplification,  we  define  the  ARMA  filter  as  a  transfer  function  operator  on  a  white  noise 
process  such  that. 


>>o[«3  = 


Biz-1) 

- r-e[n], 

acz-1) 


(5.3) 


where 


A(z  l)  —  l  +  aiz  1  +"-  +  apz~P , 

B(z  1)  =  b0+blz  1  4 - \-bQz~Q , 

and  e[n ]  is  white  noise  with  the  second  moment  being 

E{|eMf}=<r„2.  (5.6) 

From  this  definition  we  define  y^l]  as  the  output  of  the  lowpass/downsample 
procedure,  also  known  as  the  approximation  of  the  DWT.  The  goal  of  this  section  is  to 
determine  the  ARMA  filter  coefficients  of  y,[7]  as  they  relate  to  the  ARMA  coefficients 
of  the  original  signal.  Figure  V-2  shows  a  simplified  diagram  of  Method  2  where  we 
consider  a  single  stage  of  the  approximation  channel.  Once  the  procedure  to  determine 
the  ARMA  filter  coefficients  at  this  channel  is  complete,  we  will  go  back  and  consider 
the  ARMA  model  for  the  detail  channel. 


(5.4) 

(5.5) 
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Looking  again  at  Figure  V-2,  we  want  to  expand  the  innovation  process  in  the 
sense  that  we  want  to  decompose  a  signal  y0[n]  as  two  independent  signals  as  in  Figure 

V-2.  The  transfer  function  f(z_1)  is  the  lowpass  averaging  filter  described  in  equation 
(5.1)  so  that 


y,[/]  =  |j[2/]+|y[2/-i].  (5.7) 

The  sequences  e0[n]  and  ex [Z]  in  figure  V-2  are  white  noise  processes  where 
e,[Z]  is  at  half  the  sampling  rate  as  eG[n] .  Note  here  that  we  define  different  white  noise 
processes  throughout  the  derivation  and' the  terminology  can  be  confusing.  The  white 
noise  signal  ex  [Z]  is  the  final  state  of  the  multiresolution  filtering  process.  At  this  point 
there  is  no  guarantee  that  e0[ri]  and  e1  [Z]  are  uncorrelated  and  is  subject  for  future 
research.  However,  it  will  be  shown  that  the  method  does  provide  expected 
improvements  in  denoising. 


LPF 

A  /  -1  \ 
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A  (z  ) 
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_ 
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Figure  V-2  Single  Stage  of  Multiresolution  Innovation  Process 
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Starting  with  the  ARMA  model  of  the  original  signal,  we  split  the  original  signal 
into  even  and  odd  parts.  We  will  define  the  equation  now  followed  by  the  proof.  The 
goal  for  this  setup  is  to  produce  two  white  noise  processes  that  are  independent  of  one 
another.  From  the  ARMA  model  in  equation  (5.3),  we  wish  to  define  the  same  process 
decimated  into  even  and  odd  terms  giving 


CAz'1)  C0(Z-l)~ 

'  y0m  ' 

DU'1)  D(z~l) 

~ee[l] 

y0[  2/-i]_ 

z-'C'Cz-')  c,(r‘) 

join. 

D(z ")  D(z -1) 

where 


ee[l]  =  e[2l],  (5.9) 

e0[l]  =  e[2l-l],  (5.10) 

D(z~2)  =  A(z~')A{-z~x),  and  (5.11) 

Ce(z~2)  +  z~lC0(z~2)  =  B(z~1)A(-z~l) ,  (5.12) 


where 


Ce  are  the  even  coefficients  of  the  expression  B(z  ’)A(-z  ')  and 
CD  are  the  odd  coefficients  of  the  expression  B(z~l)A(-z~l) . 

By  breaking  up  the  original  signal  y0[n]  into  even  and  odd  equations,  we  get  an 

identically  equivalent  system  of  two  equations  driven  by  two  independent  white  noise 
processes  at  one-half  the  sampling  rate  of  the  original  signal.  If  in  effect  we  can 
constmct  the  sequence  y0  [ri\  from  two  independent  white  noise  processes,  then  by  the 
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properties  of  the  ARMA  model  described  in  Chapter  2,  the  inverse  processes  will  take  the 
sequence  y0[n]  and  create  two  independent  white  noise  sequences. 

The  proof  for  the  even/odd  decomposition  of  a  signal  shown  in  equation  (5.8) 
starts  with  the  ARMA  filter  determined  from  the  original  signal  y0[n]  such  that 


A(z~l ) 


(5.13) 


In  order  to  apply  the  Noble  Identity  later  on  in  the  proof,  all  transfer  functions  will 
be  described  in  terms  of  z1  instead  of  z  .  We  want  to  represent  the  ratio  of  polynomials 
in  (5.13)  by  separating  it  into  even  and  odd  powers  of  z-1  giving  us 

H(z~')  =  He(z~2)  +  Z~lH0(Z-2) .  (5.14) 

To  determine  what  He(z~l)  and  H0(z~1)  are,  we  assume  that  both  transfer 

functions  have  the  same  denominator  polynomial  as  a  function  of  z~2 .  We  accomplish 
this  by  multiplying  the  numerator  and  denominator  by  A(-z~l ) ,  i.e. 


HU'1) 


B(z~1)-  _  C(z~l) 

A(z~1)-  A(-z~l)  D(z-2)' 


(5.15) 


Recall  that  a  polynomial  can  be  decomposed  into  even  and  odd  parts  as  in 
equation  (5.14).  Performing  this  decomposition  on  the  numerator  C(z-1)  we  get  the 
polyphase  decomposition  of  H(z~l)  as 


TY  /  — 1  \  TT  /  -  2\  ,  -In  /  -2  \  C  (z.  )  -1  C0(z  ) 

H(z  )  =  He(z  )  +  z  Ha(z  )=  —  —  +z 


D(z~2) 


D{z~2) 


(5.16) 
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The  numerator  of  equation  (5.16)  can  be  thought  of  as  separating  all  the  even  and 
odd  powers  of  z-1  from  the  multiplication  of  polynomials  £(z"‘)  and  A(-z-1).  To 
show  the  denominator  is  a  polynomial  of  powers  of  z~2,  however,  takes  some 
explanation.  The  denominator  is  the  multiplication  between  polynomials  A(z-1)  and 
A(-z-1) .  We  choose  to  factor  the  polynomial  A(z_1)  in  terms  of  its  poles  p,  as 

=  (5.17) 

t=l 

where  the  total  number  of  poles  is  also  the  filter  length  P . 

In  the  same  way,  A(-z-1)  can  be  factored  as 

A(-z-‘)  =  n(l  +  p,.z-').  (5.18) 

1=1 

Substituting  equations  (5.17)  and  (5.18)  into  the  denominator  of  the  polyphase 
decomposition  equation  (5.16)  we  get 


A(z-1 )  •  A(-z_I )  =  n  (!  "  Pit'  )(1  +  PiZ~l).  (5.19) 

Z— 1 

Again,  the  single  series  product  reflects  the  fact  that  the  same  index  is  used  for 
both  expressions.  You  can  see  that  the  cross  terms  of  the  expression  inside  the  series 
product  are  zero  so  that  the  resulting  expression  yields 

A(z"‘ )  •  A(-z"‘ )  =  FI  G  “  Pt2z~2  >  =  D(z~2 )  •  (5-20) 

<=i 
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From  the  ARMA  decomposition  equation  (5.8)  we  see  that  a  signal  can  be 
synthesized  from  white  noise  processed  through  even  and  odd  parts  of  the  ARMA  model 
as  in  Figure  V-3. 


Figure  V-3  Block  Diagram  of  ARMA  Decomposition  Equation 


Now  consider  each  output  in  Figure  V-3  separately.  Starting  with  the  even  terms 
y0[2/]  of  the  original  signal  we  have  the  block  diagram  as  shown  in  Figure  V-4.  Note 

that  the  downsampling  operation  is  taken  to  each  branch  on  the  other  side  of  the 
summation.  At  this  point  we  use  the  Noble  Identity  as  in  Figure  V-5  and  we  get  the  even 
and  odd  parts  of  the  ARMA  model  in  terms  of  z-1  shown  in  Figure  V-6  [16]. 
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Figure  V-4  Block  Diagram  for  Even  Terms  of  Original  Signal 


G(z -2) 

- ► 

12 

- ►  - ► 

12 

- ► 

G(z~l) 

Figure  V-5  The  Noble  Identity 


Figure  V-6  Block  Diagram  for  Even  Terms  After  Applying  the  Noble  Identity 


From  Figure  V-6  we  see  that  the  white  noise  process  e[n]  is  split  up  into  its  even 
and  odd  terms  as  shown  in  figure  V-7. 


The  equation  to  find  the  even  terms  of  signal  y0[n]  is  then 

y0[2l]  =  He(z~')e[2l]  +  Ho(z-1)e[2l-l].  (5.21) 

Next  we  complete  the  proof  by  considering  the  odd  terms  of  the  ARMA 
decomposition  equation  (5.8)  as  in  Figure  V-8.  Following  the  same  concept  used  to 
compute  the  even  terms  of  y0[n] ,  we  apply  the  Noble  Identity  again  leading  to  the  result 

shown  in  Figure  V-9. 
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Figure  V-8  Block  Diagram  for  Odd  Terms  of  Original  Signal 


Figure  V-9  Block  Diagram  for  Odd  Terms  After  Applying  the  Noble  Identity 


Figure  V-10  Block  Diagram  to  Determine  y0  [21  - 1]  from  a  White  Noise  Process 
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Therefore,  Figure  V-10  shows  that  the  expression  for  the  odd  terms  of  signal 
;y0[n]  is  given  by 

y0[2l-l]  =  He(Z-l)e[2l  -1]+  z_1//0(z_1M2/] .  (5.22) 

Now  consider  the  even  and  odd  terms  of  the  white  noise  process.  Point  by  point, 
there  is  no  correlation  between  the  two  processes  and  therefore,  it  can  be  shown  that  a  set 
of  odd  and  even  terms  of  a  single  white  noise  process  is  completely  independent  of  each 
other.  To  reflect  this  uncorrelatedness,  we  re-define  the  odd  and  even  parts  of  a  single 
white  noise  process  as  two  completely  independent  white  noise  processes  as 

ee[l]  =  e[2l] ,  and  (5.23) 

e0[l]  =  e[2l  -1],  (5.24) 

It  follows  then  that  the  complete  model  of  the  signal  y0[n]  is  the  combination  of 

its  odd  and  even  parts  as  defined  in  equations  (5.21)  and  (5.22)  and  Figures  V-7  and  V- 
10.  Combining  these,  results  along  with  the  newly  defined  independent  white  noise 
processes  ee[l ]  and  eg[l] ,  we  get  the  system  of  equations 

y0[2Z]  =  Ht(z-l)ee[l]  +  H0(z-')e0[l]  (5.25) 

y0  [21  - 1  ]  =  He  (z_1  )e0  [/]  +  z"1  H0  (z"1  )e€  [/] ,  (5.26) 

and  a  block  diagram  described  in  Figure  V-ll.  Finally,  by  combining  equations 
(5.25)  and  (5.26)  into  matrix  form 

_y0[2Z  -1] 


ve  get 


'  He{z~x) 

~e.m 

He(z~l)_ 

‘  eom. 

(5.27) 
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which  is  the  polyphase  ARMA  of  equation  (5.8). 


Figure  V-ll  Block  Diagram  to  Model  a  Signal  from  Two  Independent  White  Noise 

Processes 
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We  have  shown  that  it  is  possible  to  model  the  even  terms  of  a  signal  separately 
from  the  odd  using  two  independent  white  noise  processes.  Based  on  this  fact,  along 
with  the  ARMA  model  given  from  the  original  signal  y0[n],  we  wish  to  determine  the 

corresponding  ARMA  model  for  yJZ] .  Recall  yJZ]  is  the  lowpass  approximation  of  the 
original  signal  yc[n]  as  described  in  the  multiresolution  operation  of  Figure  V-2. 

As  it  has  been  said  earlier,  the  multiresolution  ARMA  filter  bank  is  dependent  on 

the  filter  F(z-1).  So  from  the  averaging  filter  F(z-1)  = -^(1  + z-1),  y,[Z]  becomes 


y0m 

y0[2i-i] 


Substituting  equation  (5.27)  into  (5.28)  we  get 


y,  [/] = ^  yo  [20 + -|  y„[2/  - 1]  = 


±  1 
2  2 


*[/]  = 


"i  r 

”  He(z~l) 

~ee[l] 

_2  2_ 

_z~lH0{z-' ) 

He(z~1)_ 

e0[l]_ 

(5-28) 


(5.29) 


Carrying  out  the  multiplication  and  substituting  He(z~l )  and  H0(z  ’)  in  terms  of 

-  and  ,  -  described  in  equation  (5.16),  we  get 
D(z~l)  D{z~l ) 


C„(z-')  .C.fe-1)' 

+  Z  - 


D(z~l)  '  D(z~l) 


Cs(z2l+Cj£l 

D(z~l)  D(z~l ) 


MA- 
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Making  one  final  substitution  the  final  expression  becomes 


*P]  = 


D(z~' ) 


«.[/]  + 


D(z~l) 


e0[l]. 


where 


(5.30) 


£1 

11 

|  H- 

z~'CJz~' )]  and 

(5.31) 

c,(z-1)+c„b-1)]. 

(5.32) 

Now  consider  equation  (5.30)  as  a  difference  equation.  We  want  to  determine 
[/]  as  it  satisfies  the  equation 

N  N  M 

y,  m  -  2  /.'•«.  m +£/„”• «.  in- -  2 d,y<  p  -  mi •  (5.33) 

«=0  n=0  m=l 


where 

dp  =  coefficients  to  the  polynomial  expression  D(z_1 )  assuming  d0  =  1 , 
/„*  =  coefficients  to  the  polynomial  expression  Fe(z-1), 

//  =  coefficients  to  the  polynomial  expression  Fa(z~1 ) , 

P  ='  length  of  original  Autoregressive  (AR)  model  A{z~l ) , 

Q  =  length  of  original  Moving  Average  (MA)  model  B(z~y), 

M  =  2P  - 1 ,  and 
N  =  P  +  Q-l. 
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In  order  to  solve  this  difference  equation,  we  convert  equation  (5.30)  from  a  sum 
of  two  transfer  functions  into  a  sum  of  two  state  space  models.  Recall  the  Type  I 
controllable  canonical  state  space  form  for  the  Direct  Form  realization  of  the  transfer 
function  [13] 


H(z)  = 


b0zN+blz"-l+-bN 
zN  +  uxzN  1  ■+ +■  aN 


(5.34) 


is 


~v1[n  + 1]" 

0 

1 

v2[n  +  l] 

0 

0 

0 

0  0 

_VN  [w  +  l]_ 

1 

1 

a 

—  aN- 1 

0 

0 

1 

~a\ 


v,[n] 

'O' 

• 

v2[n] 

+ 

0 

_vA,[n]_ 

1 

x(n). 


(5.35) 


Therefore,  equation  (5.35)  is  of  the  form 

V[n  +  l]  =  AV[n]  +  Bx[n], 
and 


y[w]  —  b0aN),(bN_x  b0ciN_x),"m,(bx  i>0flj)]' 


v,[n] 

v2[n] 

vw[n] 


+  b0x[n]. 


which  leads  to 


y[n]  =  C-V[n]  +  Dx[n]. 


(5.36) 
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Since  there  are  two  transfer  functions  in  equation  (5.30),  we  apply  the  state  space 
form  twice  yielding 


xe[l  +  l]  =  Aex[l]  +  Bee[l] 

(5.37) 

yuW  =  Cexe[l]  +  Deem 

(5.38) 

x0[l  +  l)  =  A0xil]  +  B0e[l] 

(5.39) 

Vi  ,0[l]  =  C0x0[l]  +  D0e[l] 

(5.40) 

(5.41) 

where 


ce=[(feN-feodN),(rN-l-feodN_l),--,(n-feod1)] 

co=[{f°N-f%dNUf°N-x-f\dN_x),-..Xf\-f\dj\ 

De  =  feo,  and 

D0=f°  o. 
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Ae  =  Aa  because  these  matrices  are  based  on  the  same  denominator  term  D(z~l) . 
From  the  state  space  model  we  can  now  express  equation  (5.30)  in  terms  of  matrix 
variables  as 


^77 W  +  ^7777 eoW  =  {Ce (zl  -Ae)~lBe  +  De )ee [l]  +  (co (zl  -  A0 r1  B0  +  D0  )e0 [/] 


D(z~L) 


Diz'1) 


(5.42) 


Taking  the  transpose  of  expression  (5.42)  leads  to 
*[/]  =  (»/ (rf-Vr’C/  +£>,k[']  +  (s/W-A/r‘C/  +D,kPl.  (543) 

Note  that  De  and  D0  are  scalars.  There  are  clear  similarities  between  the  two 
terms  in  equation  (5.43),  which  is  advantageous  because  we  wish  to  combine  the  two 
expressions  into  one  state  space  model.  Since  Ae  and  A0  are  the  same  matrix,  the  term 
to  determine  the  eigenvalues  is  equivalent.  We  combine  both  terms  of  equation  (5.43) 
resulting  in  a  single  state  space  model  for  ^  [/] ,  giving 

x[l  +  l]  =  Ax[l]  +  Be[l]  (5.44) 

yl[l]  =  Cx[l]  +  De[l],  (5.45) 

where  the  matrices  in  the  new  state  space  model  are  defined  as 

a  =  aJ=aot, 

B  =  [c/  C/J, 

c  =  b,t=bJ, 

D  =  [De  D0  ] ,  and 
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e[l]  = 


ee[l] 

e0U]_ 


Looking  at  the  new  single  set  of  state  space  equations  (5.44)  and  (5.45),  they 
describe  a  general  stochastic  state  space  model 

x[k  + 1]  =  Ax[&]  +  v[k]  (5.46) 

z[k]  =  Cx[k]  +  \£k].  (5.47) 

An  innovation  model  can  be  determined  from  equations  (5.46)  and  (5.47)  using  a 
Kalman  Filter  as 


x[*  +  l]  =  Ai[*]  + £•?[*],  (5.49) 

z[k]=z{k]-Cx[k],  (5.50) 


where 

z[&]  is  the  observation  at  time  k, 

x[A:]  is  the  prediction  of  state  x[fc]  given  all  observations  up  to  time  k  - 1 ,  and 
z[k]  is  the  innovation,  which  is  independent  on  all  past  observations. 

The  Kalman  gain  matrix  K  and  the  covariance  matrix  P  are  computed  from  the 
Riccati  equation  where 

P  =  £{3c[k]xr[A:]} 
with  x[&]  =  x[£]  -  x[k] . 

Finally,  the  covariance  of  the  innovation  is  given  by 
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E\z[kf\=CPCT+R 


where  R  is  the  variance  of  the  white  noise  process  w[fc] . 

For  the  particular  case  of  equations  (5.44)  and  (5.45),  we  obtain 

$Ll  +  l]  =  Ax[l]  +  K-el[l],  (5.51) 

ylU]  =  Cx[l]  +  e1[l].  (5.52) 

In  order  to  solve  the  Riccati  equation,  covariance  matrices  Q ,  R  and  S  must  be 
determined.  Because  we  assume  zero-mean  white  noise,  the  covariance  matrices  are 
identically  equal  to  the  correlation  matrices.  Applying  the  general  state  space  model 
equations  (5.46)  and  (5.47)  with  the  specific  model  equations  (5.44)  and  (5.45),  we  define 
Q  as 

Q  =  Hy[&]  •  1  [*]}=  E{Be[l]eT  [l]BT  }. 

Performing  the  product  of  the  uncorrelated  2  x  N  white  noise  set,  we  note  that 


e[l]eT[l]=  e\l  [ee[l]  ea[l]]  = 


=  (J2I. 


Now  substituting  back  into  the  above  equation  and  noticing  that  all  the 
expressions  in  the  expectation  are  deterministic  we  get 


<2  =  <72BBt  . 


(5.53) 


Performing  similar  operations  to  R  and  S  we  end  up  with 
R  =  E{w[l]wr[l]}=(r2DDT , 


(5.54) 


and 

S  =  fi{v[/]  >/[/]}=  <7  W.  (5.55) 

Note  here  that  the  variance  o2  is  the  same  in  equations  (5.53)  through  (5.55) 
because  the  white  noise  processes  that  comprise  e[l]  are  independent,  yet  have  the  same 
variance.  Recall  this  variance  is  determined  from  the  initial  ARMA  model. 

For  the  final  step  in  determining  the  multiresolution  innovation  process,  consider 
the  measurement  equation  (5.52)  where  we  have  determined  all  values  in  order  to 
determine  the  innovation  ex[l]  as 


e1[l]  =  ylU]-Cm, 

where  it  can  be  shown  that  the  variance  of  ex  [/]  is 

tr12=F{ei[/]|2}=CPCr+R.  (5.56) 

It  would  be  possible  to  implement  this  multiresolution  innovation  process  in  terms 
of  a  Kalman  Filtering  recursion,  but  notice  the  state  space  model  representation  of 
equations  (5.51)  and  (5.52).  This  set  of  equations  is  also  in  Type  I  controllable  canonical 
state  space  form  for  the  Direct  Form  realization  of  a  transfer  function  [13].  This  means 
that  a  single  transfer  function  can  be  determined  from  the  state  space  model,  giving 

ym  =  [C(zl  -  A)-1  K  +  l]e,  [/] .  (5.57) 
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Finally,  it  can  be  shown  that  if  an  ARMA  model  is  minimum  phase,  then  all  its 
corresponding  multiresolution  ARMA  filters  are  minimum  phase  as  well.  From  chapter 
2,  recall  that  a  minimum  phase  filter  will  have  a  causal  stable  inverse.  From  this  fact,  we 
see  that  the  inverse  to  the  multiresolution  ARMA  model  yields  our  goal  of  obtaining  a 
multiresolution  innovation  filter  bank. 

In  a  similar  way,  we  can  determine  a  model  for  the  highpass  component,  as  in 
Figure  V-l.  Consider  now  the  ARMA  filter  coefficients  from  the  highpass  filter  channel 


where  G{z  1 )  = 


2 


-1 

2 


We  see  that  the  derivation  is  quite  similar.  In  fact,  the  only 


difference  is  the  section  where  we  consider  the  effect  ofF(z-1)  on  the  Direct  Form 
realization  of  the  transfer  function  in  equation  (5.42).  Going  back  to  equation  (5.29)  we 
substitute  the  highpass  filter  G(z-1) ,  giving 


"1  -f 

\[l] 

_2  2  _ 

z-lH0{z~l) 

ea[l]_ 

(5.58) 


Carrying  out  the  multiplication,  substituting  in  equation  (5.16),  and  grouping 
terms  yields  the  expression 


ydn  = 


GJ^) 

D(z -1) 


e.m  + 


GAz^} 

D(z~l) 


e0U], 


(5.59) 


where 


G,  <z" )  =  \[C,  (z_l )  -  z"C„(z")]  (5.60) 


and 
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(5.61) 


C/z",)  =  |[c„U_,)-C,(z-')]. 


From  this  point,  we  can  substitute  Ge(z  ‘)  and  G0(z~l)  into  equation  (5.42), 

which  is  the  state  space  form  for  the  Direct  Form  realization  of  the  transfer  function  in 
equation  (5.59). 


B.  METHOD  2  SUMMARY 

Method  1  was  a  method  to  look  at  a  whitening  filter  before  wavelet 
decomposition.  By  contrast,  Method  2  investigates  a  method  to  perform  the  whitening 
filter  after  the  multiresolution  process.  The  algorithm  can  be  divided  into  two  processes. 
The  first  one  is  the  standard  wavelet  decomposition  of  the  original  signal  similar  to  the 
DWT  on  the  whitened  signal  for  Method  1.  The  second  process  in  the  algorithm  for 
method  2  deals  with  determining  the  innovation  filter  bank.  A  set  of  ARMA  filter 
coefficients  is  necessary  for  each  wavelet  decomposition  channel,  however  all  that  is 
needed  a  priori  is  the  ARMA  model  coefficients  of  the  original  signal.  That  is  to  say,  we 
determine  the  ARMA  filter  coefficients  based  on  the  original  signal  and  from  that  model 
we  can  determine  all  follow-on  multiresolution  channels.  From  this  description,  we  can 
see  how  the  Kalman  filtering  procedure  described  in  this  chapter  is  applied  and  a  short 
description  of  the  steps  for  the  entire  process  is  shown  in  Table  V-l. 

All  the  white  noise  channels  are  uncorrelated  with  each  other  resulting  in  an 
ability  to  observe  the  innovation  process  with  varying  degrees  of  highpass  and  lowpass 
filtering.  By  way  of  analogy,  one  of  the  useful  properties  of  the  wavelet  decomposition  is 
denoising.  Since  most  useful  signal  information  occurs  in  lower  frequencies  for 
underwater  acoustic  applications,  the  wavelet  decomposition  can  leave  out  some  highpass 
information  (i.e.  details)  upon  synthesis  thereby  removing  noise  from  the  signal. 
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In  our  particular  application,  a  similar  method  of  running  a  signal  through  a 
constant-Q  filter  takes  place,  only  now  we  take  each  channel  and  process  it  through  a 
whitening  filter.  In  effect,  each  channel  has  less  and  less  high  frequency  content.  It 
would  then  stand  to  reason  that  this  multiresolution  whitening  filter  bank  might  perform 
two  desirable  functions.  First,  the  filter  bank  might  be  able  to  distinguish  low  frequency 
narrowband  signals  well  due  to  its  constant  narrowing  low  pass  filter,  providing  high 
resolution  in  low  frequencies.  Second,  since  the  channels  are  uncorrelated,  indication  of 
a  transient  signal  might  be  evident  in  one  channel  and  not  in  others,  depending  on  its 
frequency  content. 
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STEP 


DESCRIPTION 


EQUATION 


MATLAB  FUNCTION 


1 

Determine  ARMA  coefficients  of  original 
signal  and  variance  a 2 

2.3 

\b,  a\  =  cristi  _  ar(y ) 

2 

Decompose  a  filter  transfer  function  into 
even  and  odd  parts  for  Low  Pass  Filter 

F(z~ *) 

5.15,  5.16 

\Ce,C0,D]=-  polyd ( num , den ) 

3 

Determine  Fe(z  *)  and  F0(z  ') 

5.31,  5.32 

l^n+l  ’  ^n+l  ’  an+ 1  ’  fln+l .  ~ 

model2{bn,an ) 

4 

Determine  the  Type  I  Controllable 

Canonical  Form  for  He(z~x)  and#0(z_1) 

5.35,  5.36 

tf2ss(F„D) 

5 

Combine  even  and  odd  processes  to  one 
state  space  model 

5.44,  5.45 

l^n+1  ’  b„+ 1  )  ^„+l  J  ^n+l  j  — 
model2(bn,an) 

6 

Determine  Q,  R  and  S 

5.53  -  5.55 

Q  =  <72BBt  ,  R  =  cr2DDT 

S  =  o2BDt 

7 

Determine  Kalman  Gain  and  State 
Covariance  Matrix  from  Riccati  Equation 

[14] 

[K,P]  = 

mydlqe(A,  G,C,Q,R,S ) 

8 

Determine  new  variance  a 2 

5.56 

a2  =CPCT  +R 

9 

Determine  next  multiresolution  filter 
transfer  function 

5.57 

[num,  den]  = 
ss2tf(A,K,C,l) 

Bfl 

Repeat  steps  2  -  12  for  HPF  G(z-1) 

5.60,  5.61 

Included  in  Step  5 

H 

Perform  2  -  13  to  the  desired  number  of 
stages 

meth_0102 

12 

Run  the  signal  through  the  multiresolution 
whitening  filter  bank  as  in  Figure  V-l 

Figure  V-l 

meth_0102 

Table  V-l  Process  to  Determine  the  Multiresolution  Whitening  Filter  Bank 
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VL  RESULTS 


Now  that  we  have  set  a  framework  for  two  methods  of  multiresolution  filtering,  it 
is  necessary  to  compare  the  techniques  against  the  benchmark  signal  described  in  Chapter 
1.  Method  1  can  be  described  as  pre-whitening  a  signal  before  multiresolution  analysis, 

t 

which  is  the  Discrete  Haar  Wavelet  decomposition.  By  contrast,  a  signal  processed 
through  Method  2  uses  the  same  two  processes  in  opposite  order.  Again,  the  benchmark 
signal  is  a  fan  blowing,  which  is  representative  of  the  colored  noise  environment  of  the 
heavily  -  trafficked  coastal  waters  because  it  contains  more  of  the  lower  frequencies  due 
to  a  rotating  blade  through  a  fluid. 

The  test  signal  used  in  this  section  is  the  first  three  seconds  of  the  benchmark 
signal  unless  otherwise  specified.  The  assumption  is  there  are  no  other  transients  in  this 
regime  except  for  the  synthetic  signal  that  starts  after  the  first  second.  A  wide  range  of 
synthetic  transients  is  tested  from  narrowband  single  frequency  sinusoids  to  broadband 
white  noise.  The  object  of  this  section  is  to  compare  established  detection  techniques 
with  Methods  1  and  2  using  a  wide  range  of  transients. 

Since  the  STFT  is  the  current  standard  for  underwater  acoustic  analysis,  all  results 
will  begin  with  a  spectrogram  followed  by  the  whitening  filter  process.  In  this  way,  the 
viewer  can  receive  an  intuitive  “feel”  for  comparing  methods.  Introducing  the  new 
detection  methods,  as  we  shall  see,  results  in  a  higher  confidence  in  positively  identifying 
the  transient. 

Finally,  for  ease  of  comparison.  Method  1  and  2  results  will  be  presented 
concurrently  at  each  of  the  three  transient  tests.  Because  of  the  multiresolution  process, 
each  method  has  multiple  channels  to  display.  Since  not  all  channels  will  indicate  the 
presence  of  the  signal,  all  figures  left  out  are  assumed  to  be  clutter.  For  convenience,  a 
small  diagram  of  each  process  will  precede  every  figure  with  its  respective  filter  path 
highlighted.  Knowing  which  channel  is  being  displayed  is  important  when  analyzing  the 
effect  of  the  filtering  process  on  the  transient  signal. 
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A. 


BENCHMARK  TEST 


The  method  comparison  will  test  three  different  transient  types.  The  first  two  are 
narrowband  transients  at  opposite  sides  of  the  frequency  spectrum  while  the  last  test  will 
be  the  more  general  broadband  type.  We  will  start  with  examining  the  narrowband 
transient  signals  and  make  adjustments  to  Methods  1  and  2  as  necessary. 

1.  Narrowband  Transient:  High  Frequency  Sinusoid 

Note  that  the  maximum  testable  frequency  is  half  the  sampling  frequency  (just 
over  5  kHz)  as  the  benchmark  signal  is  sampled  at  11.025  kHz.  Here,  we  arbitrarily 
choose  a  5  kHz  sinusoidal  transient  with  a  duration  of  .1  second.  The  spectrogram  in 
Figure  VI- 1  shows  the  transient  is  emerging  from  the  background  but  is  barely 
recognizable  audibly  because  the  amplitude  is  small.  From  this  signal  representation  it 
can  be  seen  that  since  the  benchmark  contains  little  frequency  in  the  region  above  2.5 
kHz,  a  transient  signal  above  this  frequency  can  be  easily  recognizable,  assuming  the 
amplitude  is  high  enough. 

Figure  VI-2  shows  a  simple  diagram  of  Method  1  in  order  to  keep  track  of  the 
multiple  outputs  in  the  analysis  tree  followed  by  Figure  VI-3,  which  shows  the  innovation 
process  on  the  benchmark  signal.  Again,  in  this  representation  the  transient  is  noticeable 
between  the  vertical  dashed  lines.  Notice,  though,  the  mean  and  variance  of  the  noise 
floor  compared  to  Figure  VI-5,  which  is  the  first  stage  of  the  multiresolution  process.  In 
this  procedure,  the  whitened  signal  is  now  processed  through  a  lowpass  filter  and  a 
highpass  filter,  which  are  the  Haar  Wavelet  coefficients. 

Looking  at  the  results  of  the  first  stage  of  the  Haar  Wavelet  decomposition  shown 
in  Figure  VI-5,  the  signal  through  the  lowpass  filter  provides  no  useful  information,  but 
the  highpass  filter  shows  the  transient  very  clear  and  with  less  variance  on  the  noise  floor. 

As  a  minor  digression,  it  is  important  to  recall  that  the  signal  passed  through  the 
Haar  Wavelet  filters  will  result  in  the  Haar  Wavelet  coefficients,  which  satisfy  Parseval’s 
theorem  similar  to  the  properties  of  the  coefficients  determined  from  the  Fourier 
Transform.  In  this  application,  the  actual  numerical  power  given  for  each  scale  is  not 
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important,  only  their  relative  values  and  so  all  the  signals  can  be  normalized.  For  this 
application,  the  fundamental  properties  of  the  Discrete  Wavelet  Transform  are  only 
important  in  terms  of  defining  the  filtering  process.  To  effectively  analyze  the  results,  we 
simply  view  the  multiresolution  filtering  process  in  terms  of  how  the  filters  effect  the 
signal  in  the  frequency  domain.  What  we  have  shown  in  Figures  VI-3  and  VI-5  is  a 
denoising  process,  which  is  one  of  the  many  uses  of  wavelet  transforms. 
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Figure  VI-2  Method  1  Analysis  Tree  (w[n]) 


Figure  VI-3  Transient  1:  Whitening  Filter 


Figure  VI-4  Method  I  Wavelet  Analysis  Tree  (  Ax  and  Dx ) 


LOWPASS  FILTER  STAGE  1 


Figure  VI-5  Transient  1:  Method  1,  Stage  1  (A]  and  Dx) 


For  the  high  frequency  sinusoid,  the  transient  will  not  appear  in  any  follow-on 
stages  of  the  lowpass  filtering  processes  (approximations).  At  this  point,  we  are  left  with 
a  successfully  detected  signal  in  the  high  pass  channel  of  the  first  stage.  For  this 
scenario,  we  can  see  the  transient  above  the  noise  floor,  but  recall  this  signal  can  also  be 
detected  from  the  spectrogram.  Now,  suppose  we  can  reduce  the  amplitude  of  the 
transient  signal  to  the  point  where  it  cannot  be  easily  identified  in  either  the  spectrogram 
or  the  whitening  filter.  Can  the  denoising  process  of  Method  1  provide  better  results? 
The  answer  is  shown  in  the  next  set  of  figures. 

The  key  to  successful  denoising  on  the  high  frequency  end  of  the  spectrum  is  to 
filter  the  details  (high  pass  channel)  in  the  same  manner  as  the  approximations  (low  pass 
channel).  Knowing  that  we  can  remove  the  signal  from  the  noise  on  the  low  pass 
channel,  we  apply  the  same  wavelet  decomposition  at  every  output.  This  description  is 
known  as  the  standard  Wavelet  Packet  framework,  which  has  been  applied  to  a  wide 
range  of  applications  such  as  denoising  and  data  compression.  Wavelet  packet 
processing  allows  for  a  more  complex  and  flexible  analysis  because  both  the  details 
(highpass  channel)  and  the  approximations  (lowpass  channel)  are  decomposed,  as  shown 
in  Figure  VI-7.  In  terms  of  wavelet  analysis,  the  wavelet  packet  framework  allows  for  a 
wider  range  of  bases  from  which  you  can  choose  the  best  representation  [15]. 

For  our  application,  the  flexibility  of  the  wavelet  packet  allows  for  greater 
narrowband  detection  sensitivity  in  both  high  and  low  frequencies.  Take  the  high 
frequency  transient  for  example.  Let’s  reduce  the  amplitude  of  the  5  kHz  sinusoid  from 
.005  to  .003.  To  quantify  this  reduction,  audibly  detecting  this  transient  signal  is  like 
straining  your  ears  to  listen  for  the  softest  high  frequency  sounds  during  an  audiogram 
test.  On  top  of  that,  we  still  have  the  noise  of  the  fan  blowing.  Figure  VI-6  shows  the 
spectrogram  of  the  benchmark  signal  with  the  barely  noticeable  transient  at  5  kHz,  which 
could  be  classified  as  noise.  Now  we  add  the  high  pass  channel  decomposition  and  find 
that  several  channels  have  successfully  reduced  the  noise  floor  so  that  the  signal  can  be 
detected.  Because  of  the  branching  effect  of  the  wavelet  packet,  a  three-stage  framework 
will  yield  15  plots.  In  order  to  reduce  the  number  of  plots,  only  the  graphs  that  detect  the 
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signal  are  plotted.  Figure  VI-8  shows  the  whitening  filter  representation  of  the 
benchmark  signal.  Notice  that  the  signal  is  completely  imbedded  within  the  noise  such 
that  it  cannot  be  identified.  Performing  the  first  stage  of  the  wavelet  decomposition  as 
before  yields  a  slight  indication  of  the  transient  (Figure  VI-10),  but  now  we  perform  the 
filtering  process  on  the  details  and  notice  the  signal  clearly  emerges  from  the  noise  floor 
(Figure  VI-12).  Performing  the  filtering  process  one  more  time  yields  an  even  better 
result  (Figure  VI- 14).  Realistically,  the  transient  might  have  been  positively  identified  on 
the  spectrogram.  However  with  the  wavelet  denoising  process,  the  added  signal  detection 
capabilities  might  add  to  the  measure  of  confidence. 

At  first  glance,  it  would  appear  that  the  5  kHz  signal  should  be  detected  in  the  bin 
to  the  far  right  marked  DDD3  because  the  bandpass  filter  in  this  region  contains  the 
frequency  of  the  transient.  Note  in  Figures  VI-11  to  VI-14  that  the  signal  appears  in  the 
bandpass  filter  region  AD2  and  AADZ .  The  reason  for  this  could  be  due  to  aliasing  that 

comes  from  the  downsampling  process  is  subject  for  further  research.  Also,  since  the 
highpass  and  lowpass  filters  of  the  Haar  wavelet  do  not  have  a  sharp  dropoff  at  the  cutoff 
frequency,  it  is  possible  to  have  part  of  a  high  frequency  transient  signal  to  appear  in  the 
low  frequency  channel  and  vice  versa.  Both  these  issues  are  subject  for  further  research. 
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FREQUENCY  (Hz) 


Figure  VI-6  Spectrogram  of  Modified  Transient  Is  Reduced  Amplitude 
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Figure  VI-7  Wavelet  Packet  Framework  (w[n]) 


Figure  VI-8  Modified  Transient  1:  Whitening  Filter 


85 


Figure  VI-9  Wavelet  Packet  Framework  ( Ax  and  D{ ) 
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Figure  VI-12  Modified  Transient  1:  Method  1,  Stage  2  ( AD2  and  DD2) 
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Figure  VI-13  Wavelet  Packet  Framework  ( AAD^  and  DAD 3 ) 


Figure  VI-14  Modified  Transient  1:  Method  1,  Stage  3  ( AAD3  and  DADi ) 
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Turning  now  to  Method  2,  recall  that  the  input  signal  is  whitened  after  the 
multiresolution  Wavelet  decomposition  as  described  in  Figure  V-l.  Throughout  the 
derivation,  we  used  a  slightly  modified  filter  bank  from  the  Haar  Wavelet  coefficients  as 
described  in  equation  (5.7).  As  mentioned  before,  the  only  difference  is  a  constant  factor 
and  since  the  results  are  normalized,  the  shape  of  the  signal  is  exactly  the  same.  For  the 
sake  of  consistency,  the  actual  Haar  Wavelet  filter  coefficients  were  used. 

Unfortunately,  Method  2  did  not  detect  the  narrowband  high  frequency  transient 
signal  at  amplitude  .005.  Looking  at  the  lowpass  filter,  notice  that  it  averages  two 

adjacent  data  points  and  weights  them  by  a  factor  of  V2 .  For  the  case  of  the  5kHz 
signal,  perhaps  the  transient  cannot  be  seen  in  the  lowpass  channel  because  of  its  high 
frequency  content  and  cannot  be  seen  in  the  highpass  channel  because  the  transient  is 
averaged  in  with  the  noise.  Nonetheless,  to  show  that  the  method  is  still  capable  of 
detection,  it  can  be  shown  experimentally  that  if  the  amplitude  of  the  transient  signal  is 
increased  to  .01,  the  transient  signal  emerges  from  the  noise. 

Although  the  multiresolution  innovation  process  did  not  perform  as  well  as  the 
wavelet  packet  denoising  process  at  5  kHz,  it  is  still  a  viable  detection  scheme.  Indeed, 
we  shall  see  other  examples  where  Method  2  provides  the  best  transient  detection 
capabilities. 
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2.  Narrowband  Transient:  Low  Frequency  Sinusoid 

At  high  frequencies,  the  spectrogram  works  quite  well  in  detecting  narrowband 
transients  because  the  stationary  colored  noise  model  has  very  little  high  frequency 
content.  Therefore  any  added  frequency  content  is  easily  recognizable.  Now  consider  a 
transient  signal  buried  in  the  low  frequency  band  between  0  and  300  Hz.  Since  the 
colored  noise  signal  has  its  frequency  content  in  this  range,  a  signal  cannot  be 
distinguished  from  the  background  noise.  Since  both  Methods  1  and  2  are  based  on  a 
signal’s  statistics  in  the  time  domain  and  not  on  the  Fourier  Transform  in  the  frequency 
domain,  they  should  detect  low  frequency  transients  equally  as  well.  At  each  stage  of  the 
wavelet  decomposition,  the  low  frequency  channel  continues  to  be  passed  through  a 
lowpass  filter,  effectively  cutting  the  frequency  content  of  the  remaining  signal  in  half. 
Continuing  this  process  provides  infinite  resolution  in  the  low  frequency  channels,  which 
again  is  one  of  the  nice  properties  of  the  Discrete  Wavelet  Transform.  For  this  reason, 
Method  2  should  perform  at  least  as  well  as  Method  1  in  the  high  frequency  narrowband 
case. 

We  choose  an  arbitrary  sinusoidal  signal  at  50  Hz  and  conduct  a  similar  analysis. 
Figure  VI- 15  shows  the  spectrogram  of  the  benchmark  signal  along  with  the  imbedded  50 
Hz  transient  signal.  It  is  hard  to  tell  whether  or  not  the  transient  can  be  detected  in  this 
region  because  it  could  easily  be  identified  as  a  spurious  point.  The  transient  clearly 
cannot  be  detected  with  just  the  whitening  filter  (Figure  VI- 17)  and  not  even  audibly. 
Starting  with  Method  1,  we  see  the  signal  emerging  on  the  low  pass  channel  in  Figure  VI- 
19.  Because  the  frequency  is  so  low,  the  signal  becomes  more  pronounced  with  each 
successive  lowpass  filter  as  shown  in  Figures  VI-20  through  23. 
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Figure  VI-15  Spectrogram  of  Transient  2:  50  Hz  Sinusoid 


91 


Figure  VI-16  Wavelet  Packet  Framework  (w[n]) 


Figure  VI-17  Transient  2:  Whitening  Filter 
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Figure  VI-18  Wavelet  Packet  Framework  (A,  and  £>,) 


LOWPASS  FILTER  (APPROXIMATION) 


Figure  VI-19  Transient  2:  Method  1,  Stage  1  (A!  and  D.) 


Figure  VI-20  Wavelet  Packet  Framework  (AAj  and  DA2) 


LOWPASS  FILTER  (APPROXIMATION) 


Figure  VI-21  Transient  2:  Method  1,  Stage  2  ( AA2  and  DA2) 


Figure  VI-22  Wavelet  Packet  Framework  ( AAA3  and  DAA3 ) 


Figure  VI-23  Transient  2:  Method  1,  Stage  3  ( AAA 3  and  DAA3 ) 
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Method  2  results  for  the  low  frequency  signal  are  just  as  revealing.  Unlike  the 
Wavelet  Transform  on  the  white  noise  process,  the  colored  noise  signal  does  not  have 
equal  frequency  content  at  all  frequencies  (except  for  the  transient)  and  therefore  the 
narrowband  transient  will  not  fall  on  one  side  of  the  highpass/lowpass  filtering  process. 
Indeed,  the  signal  may  show  up  in  as  many  as  all  of  the  channels  as  is  the  case  here, 
which  in  and  of  itself  might  provide  added  redundancy  since  the  channels  are 
independent  of  each  other.  For  the  lowpass-filtering  channel,  the  benchmark  signal  is 
downsampled  and  averaged.  As  it  averages  adjacent  points,  the  process  essentially 
reduces  the  noise  level  while  maintaining  a  strong  indication  of  the  presence  of  a  signal. 
For  the  highpass  filter,  the  difference  between  adjacent  points  is  computed  and  the 
combination  of  the  two  processes  form  a  basis  set  for  the  original  signal.  Again,  this 
highpass  process  might  prove  valuable  since  it  is  independent  of  the  lowpass  process. 
Figure  VI-25  shows  the  signal  clearly  emerging  in  both  channels.  With  each  successive 
multiresolution  innovation  process,  the  noise  floor  is  reduced  and  the  signal  becomes 
more  pronounced  (Figures  VI-26  to  29). 
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Figure  VI-25  Transient  2:  Method  2,  Stage  1(AX  and  Dx) 
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Figure  VI-26  Multiresolution  Innovation  Stage  2  {AA2  and  DA2) 


Figure  VI-27  Transient  2:  Method  2,  Stage  2  ( AA2  and  DA2 ) 
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Figure  VI-28  Multiresolution  Innovation  Stage  3  ( AAA3  and  DAA.) 


LOWPASS  FILTER  (APPROXIMATION) 


Figure  VI-29  Transient  2:  Method  2,  Stage  3  ( AAA3  and  DAA3) 
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3.  Broadband  Transient:  Elvis 

We  have  studied  the  effects  of  Methods  1  and  2  on  narrowband  transients.  In  the 
case  of  underwater  acoustics,  it  is  quite  possible  that  the  signal  we  wish  to  detect  is  a 
single  frequency  like  an  active  sonar  pulse  or  a  unique  single  frequency  tonal.  On  the 
other  hand,  a  broadband  signal  such  as  the  splash  of  a  mine  hitting  the  water  or  the  clang 
of  a  metallic  object  is  also  possible. 

From  the  narrowband  transient,  we  now  consider  a  more  general  case  where  the 
expected  signal  has  a  wide  frequency  band.  For  this  test,  the  frequency  content  must  be 
completely  arbitrary,  so  for  no  apparent  scientific  reason  the  arbitrary  colored  noise 
signal  will  be  Elvis  singing  that  great  Christmas  classic,  “Blue  Christmas”.  Besides,  it  is 
just  as  likely  to  find  Elvis  in  the  deep  blue  ocean  as  it  is  to  find  him  anywhere  else. 

In  order  to  see  the  frequency  content  of  the  narrowband  transient  signal,  a 
spectrogram  of  the  transient  alone  is  provided  in  Figure  VI-30(a).  Notice  that  Elvis  can 
really  hit  those  low  notes.  The  strength  of  the  frequency  content  is  large  in  the  range 
between  0  and  1  kHz,  then  tapers  off  as  it  approaches  5  kHz.  Adding  Elvis  to  the  colored 
noise  model  as  in  Figure  VI-30(b),  notice  that  the  signal  can  barely  be  seen,  and  that  only 
because  we  have  other  spectrograms  of  the  same  data  with  which  to  compare.  Without 
the  redundant  benchmark  data,  “Blue  Christmas”  is  reduced  to  mere  noise.  We  perform 
the  innovation  process  on  the  benchmark  signal  in  Figure  VI-32  and,  no  doubt  about  it, 
Elvis  is  clearly  alive.  Performing  the  denoising  on  the  innovation  as  in  Method  1 
confirms  this  to  be  true,  however  the  signal  is  less  pronounced.  Also  notice  that  the 
signal  shows  up  -  in  part  -  on  both  highpass  and  lowpass  channels  at  the  same  stage. 
This  indicates  that  the  broad  range  of  frequency  content  of  Elvis’  golden  voice  is  divided 
into  highpass  and  lowpass  regions.  The  problem  with  Method  1  here  is  that  along  with 
the  reduction  of  noise  with  each  filter  stage  comes  the  reduction  of  the  signal  power  as 
well  depending  on  its  frequency  content.  In  other  words,  the  success  of  Method  1  is  now 
dependent  on  frequency  content. 

The  signal  had  the  same  effect  when  processed  through  Method  2,  only  worse. 
Only  the  first  stage  provided  any  indication  of  Elvis’  location  as  in  Figure  VI-44.  Here 
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again,  the  broadband  signal  could  have  been  averaged  in  with  the  noise  so  that  the 
innovation  process  could  not  distinguish  Elvis  from  a  fan. 
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Figure  VI-30  Spectrogram  of  Elvis  Without  Benchmark  Colored  Noise 
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Figure  VI-31  Spectrogram  of  Transient  3:  Elvis 
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Figure  VI-32  Wavelet  Packet  Framework  (w[n]> 


Figure  VI-33  Elvis  Through  a  Whitening  Filter 
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Figure  VI-34  Wavelet  Packet  Framework  (Al  and  Dx) 


Figure  VI-35  In  Search  of  Elvis:  Method  1,  Stage  1  (At  and  Dx) 
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Figure  VI-37  In  Search  of  Elvis:  Method  1,  Stage  2  ( AA^  arid  DA2) 
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Figure  VI-38  Wavelet  Packet  Framework  (ADA3  and  DDA3 ) 
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Figure  VI-40  Wavelet  Packet  Framework  ( AD2  and  DD2) 
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Figure  VI-41  In  Search  of  Elvis:  Method  1,  Stage  3  ( AD2  and  DD2) 


107 


Figure  VI-42  Wavelet  Packet  Framework  (ADD3  and  DDD3 ) 
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Figure  VI-43  In  Search  of  Elvis:  Method  1,  Stage  3  ( ADD3  and  DDD3 ) 


Figure  VI-44  Multiresolution  Innovation  Stage  1(A1  and  Dx) 


Figure  VI-45  In  Search  of  Elvis:  Method  2,  Stage  l(A1  and  Dx) 
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B.  SUMMARY  OF  OBSERVATIONS 


We  have  presented  a  representative  set  of  transient  signals  tested  on  two  detection 
schemes.  This  section  aims  to  provide  a  broad  brushstroke  of  the  effectiveness  of  the 
Methods  as  they  are  compared  to  the  spectrogram  and  the  whitening  filter.  To 
accomplish  this  comparison,  a  measure  of  successfulness  must  be  defined.  Since  the 
location  and  duration  of  the  transient  signal  is  known,  an  estimate  of  the  signal  power  and 
the  noise  power  can  be  calculated  as  the  Signal-to-Noise  ratio.  Since  our  signal  model 
has  additive  noise,  the  signal  power  is  simply  the  power  above  the  noise  threshold 
divided  by  the  noise  power  where  both  power  terms  can  be  calculated  from  the  whitened 
signal.  For  Methods  1  and  2,  whose  multiple  channels  will  not  always  provide  the  best 
transient  detection,  only  the  channel  with  the  maximum  SNR  will  be  taken.  The 
assumption  here  is  that  there  is  a  way  of  always  determining  the  best  channel  within  a 
given  method  scheme. 

For  a  signal  whose  power  is  twice  the  noise  power,  the  SNR  will  be  3  dB.  This 
value  will  be  our  threshold  detection  mark  where  any  SNR  above  this  value  is  assumed  to 
have  positively  detected  the  transient.  Any  value  below  this  mark  is  assumed  to  not  be 
able  to  identify  the  transient  even  though  a  visual  representation  might  indicate 
otherwise. 

As  in  all  experiments,  the  SNR  measure  is  only  an  approximation.  In  our  case, 
the  same  exact  transient  signals  within  the  colored  noise  model  is  used  to  compute  the 
SNR  and  so  the  results  among  the  different  methods  may  be  compared  to  one  another.  It 
is  fair  to  assume,  however,  that  the  results  from  this  experiment  are  indicative  of  many 
types  of  transient  signals  hidden  within  an  arbitrary  colored  noise  scheme.  The  first 
experiment  deals  with  broadband  signal  detection  followed  by  a  test  on  narrowband 
signal  detection. 
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1.  Benchmark  Comparison 

The  first  experiment  begins  with  Elvis.  Taking  the  signal,  we  increase  the 
normalized  volume  from  0  to  1  and  determine  its  SNR  from  the  whitening  filter  alone, 
Method  1  and  Method  2.  Figure  VI-45  shows  the  results.  Here,  it  can  be  seen  that 
Method  1  reaches  the  3-dB  threshold  first,  followed  closely  by  Method  2  and  the 
benchmark  whitening  filter.  As  Elvis  gets  louder,  though.  Method  1  reaches  a  higher 
maximum  SNR  value  of  25  dB  while  the  other  two  methods  are  just  about  identical  as 
they  reach  a  steady  state  value  of  about  20  dB.  The  results  of  the  experiments  show  that 
the  noise  level  does  indeed  reduce  faster  than  the  signal  and  so  a  better  detection  method 
was  found  even  for  arbitrary  broadband  signals.  On  the  other  hand,  comparing  any  of 
these  results  with  the  threshold  reveals  that  Method  1  might  reach  detection  status  before 
the  other  methods,  but  not  by  much.  What’s  more,  the  spectrogram  might  be  able  to  pick 
up  the  signal  before  any  of  the  three  signals,  but  again  there  must  be  some  frequency 
content  outside  of  the  frequencies  defined  by  the  colored  noise. 

Another  broadband  test  was  conducted  where  the  frequency  content  expanded  the 
entire  spectrum.  A  similar  test  was  conducted  where  the  broadband  signal  was  white 
noise  and  the  results  are  found  on  Figure  VI-46.  Here  again,  Method  1  broke  out  with  the 
highest  maximum  SNR  at  40  dB.  Since  all  three  methods  reached  the  3  dB  threshold  at 
about  the  same  amplitude,  all  three  methods  perform  about  the  same.  Also,  we  are 
guaranteed  to  have  the  spectrogram  detect  the  transient  at  some  point  close  to  the  other 
methods  and  so  any  of  these  processes  are  equally  effective. 
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Figure  VI-46  SNR  Comparison  on  Colored  Noise  Transient 


Figure  VI-47  SNR  Comparison  on  White  Noise  Transient 
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The  final  comparison  is  with  the  narrowband  transient.  For  this  test  we  vary  both 
frequency  and  amplitude  of  a  standard  .  1  second  transient  and  process  this  signal  through 
the  whitening  filter.  Method  1  and  Method  2.  The  resulting  3-dimensional  image  reveals 
some  interesting  peaks  and  valleys.  Figure  VI-47  shows  the  result  from  the  whitening 
filter.  Comparing  this  plot  with  the  results  from  Method  1  (Figure  VI-48)  we  see  right 
away  that  it  is  the  same  basic  shape,  only  the  SNR  has  increased  by  about  the  same 
amount  over  the  entire  frequency-amplitude  plane.  Now  recall  that  the  threshold  is  3  dB 
and  notice  how  well  the  wavelet  denoising  process  worked.  In  the  same  manner,  Method 
2  produced  very  good  results  in  the  low  frequency  region. 

We  wish  to  further  quantify  the  results  of  the  comparison  test.  Taking  the  3- 
dimensional  plot,  we  slice  the  graph  along  the  frequency-amplitude  plane  at  a  constant 
level  of  3  dB.  Looking  down  on  the  slice,  we  can  see  clearly  the  regions  where  the  SNR 
is  above  3  dB,  which  is  the  detection  region.  Figure  VI-49  shows  the  modest  narrowband 
transient  region  for  the  whitening  process.  With  the  exception  of  a  small  strip,  the  lower 
frequency  SNRs  remain  below  the  threshold,  especially  between  0  and  300  Hz  where 
there  is  no  SNR  above  3dB.  Now  Comparing  Method  1  to  the  whitening  filter  detection 
region  (Figure  VI-50)  we  notice  right  away  the  larger  areas  in  which  the  SNR  is  above 
the  threshold.  To  quantify,  Table  VI- 1  shows  that  only  about  30%  of  the  entire 
frequency-amplitude  plane  is  above  the  threshold.  With  the  Haar  Wavelet  denoising 
method,  the  percent  detection  over  the  entire  region  increased  to  74%.  In  the  low 
frequency  region  between  0  and  2  kHz,  the  detection  region  increased  from  29.7%  to 
78.3%.  The  results  from  Method  2  as  seen  in  Figure  VI-51  show  an  equally  strong 
performance  in  the  low  frequency  regions  and  as  mentioned  earlier,  it  has  limitations  in 
the  high  frequency  regime. 

These  results  seem  very  good  compared  to  the  whitening  filter,  now  compare  the 
methods  against  the  spectrogram.  The  spectrogram  can  detect  most  signals  above  an 
amplitude  of  .1  and  frequencies  above  300  Hz.  Recall  that  the  whitening  filter  did  not 
have  any  of  its  detection  region  below  300  Hz  so  with  these  two  methods  combined,  low 
frequency  transients  would  go  unnoticed.  Looking  back  at  Figures  VI-50,  VI-51  and 
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Table  VI-1,  we  notice  that  the  detection  region  above  the  3  dB  mark  is  now  almost  50% 
of  the  total  area  in  this  region. 
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PERCENT  DETECTION  IN  FREQUENCY  RANGE 


FILTERING 

TECHNIQUE 

0  -  5  (kHz) 

0  -  2  (kHz) 

0  -  300  (Hz) 

WHITENING  FILTER 

30.3 

30.3 

0.0 

METHOD  1:  HAAR 

WAVELET  DE¬ 
NOTING 

73.9 

79.6 

50.5 

METHOD  2: 

MULTIRESOLUTION 

INNOVATIONS 

66.1 

81.1 

51.8 

Table  Vl-1  Narrowband  Transient  Detection  Capability  Comparison  Over  a 

Threshold  of  3  dB 
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Figure  VI-49  Narrowband  Transient  SNR  for  Haar  Wavelet  Denoising 
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TRANSIENT  AMPLITUDE 


SNR  (dB) 


TRANSIENT  FREQUENCY  (Hz) 


Figure  VI-50  Narrowband  Transient  Detection  Region  For  Whitening  Filter 

(Threshold  =  3  dB) 
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Figure  VI-51  Narrowband  Transient  Detection  Region  For  Haar  Wavelet  Denoising 

(Threshold  =  3  dB) 
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Figure  VI-52  Narrowband  Transient  Detection  Region  For  Whitening  Filter 

(Threshold  =  3  dB) 
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Figure  VI-53  Narrowband  Transient  Detection  Region  For  Multiresolution 
Innovation  Process  (Threshold  =  3  dB) 
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2.  Strengths 

The  most  useful  property  of  the  wavelet  denoising  process  and  the  multiresolution 
innovation  process  are  their  abilities  to  detect  narrowband  signals  in  the  low  frequency 
regime,  particularly  in  between  0  and  300  Hz.  This  property  is  highly  desirable  since  it 
has  been  shown  that  neither  the  spectrogram  nor  the  whitening  filter  is  effective  in  this 
frequency  region. 

The  Haar  Wavelet  denoising  routine  works  well  at  high  frequencies  too  because 
of  its  wavelet  packet  framework.  Here,  the  high  frequency  details  are  filtered  too, 
creating  greater  flexibility  by  isolating  more  frequencies. 

Since  the  multiple-channel  outputs  of  both  Method  1  and  Method  2  are 
independent  of  each  other,  it  is  likely  that  more  than  one  channel  will  indicate  the 
presence  of  a  transient  signal.  This  property  provides  a  measure  of  redundancy  in 
positively  identifying  a  transient  signal. 

Both  methods  are  simple  to  compute  numerically;  the  DWT  is  computed  by  a 
simple  filtering  process  as  the  number  of  computations  is  comparable  to  current  methods 
that  use  the  FFT.  Since  the  size  of  the  data  set  is  reduced  by  a  factor  of  2  at  each  stage, 
the  number  of  computations  is  reduced  by  the  same  amount.  This  property  is  desirable 
because  we  know  that  any  follow-on  stages  will  incur  less  and  less  computations. 

Both  methods  lend  themselves  to  an  efficient  modular  processing  routine  with  a 
user-friendly  GUI  display.  Take  for  example  the  wavelet  packet  framework.  The  user 
might  designate  the  number  of  stages  to  explore  and  at  any  time,  can  add  another  stage 
without  having  to  go  back  and  compute  the  entire  wavelet  tree.  The  same  is  true  for  the 
multiresolution  whitening  process  because,  like  the  wavelet  decomposition,  every 
successive  ARMA  filter  set  is  only  based  on  the  previous  stage  and  does  not  need  to  be 
computed  from  the  very  first  stage  again. 
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3. 


Weaknesses 


Both  wavelet  denoising  and  multiresolution  innovation  procedures  attempt  to 
reduce  the  noise  variance  by  a  series  of  filtering  processes.  For  a  signal  that  is  isolated  to 
a  small  range  of  frequencies,  both  methods  will  perform  quite  well  since  the  filtering 
process  has  the  ability  to  cut  off  a  sizeable  amount  of  extraneous  frequency  content.  By 
contrast,  a  broadband  signal  will  have  part  of  its  frequency  content  cut  off  at  each 
filtering  stage.  The  effect  will  be  a  reduction  of  signal  power  along  with  a  reduction  in 
noise  power.  Even  with  the  signal  power  reduction,  the  results  of  the  previous  section 
show  that  Methods  1  and  2  were  able  to  reduce  the  noise  floor  as  the  amplitude  of  the 
transient  signal  increased.  Comparing  these  methods  to  the  spectrogram,  however,  we 
see  that  the  transient  will  appear  as  a  vertical  line  at  amplitude  less  than  .05,  which  is 
about  when  the  signal  emerges  from  the  noise  in  the  multiresolution  filtering  process. 
The  advantage  of  the  spectrogram  for  broadband  signals  is  its  color  representation,  which 
essentially  adds  an  extra  dimension  to  detecting  the  transient.  The  only  exception  to  this 
observation  is  when  the  transient  contains  all  or  part  of  the  same  frequencies  as  the 
background  noise,  in  which  case  the  signal  is  camouflaged  in  the  spectrogram. 

Inherent  in  the  wavelet  decomposition  process  is  a  reduction  of  resolution  in  the 
time  domain.  As  described  in  Chapter  3,  to  perform  a  DWT  on  a  signal  it  is 
downsampled  as  part  of  the  process.  By  this  action,  the  signal  in  the  time  domain  loses 
some  detail  at  each  stage.  Figure  VI-53  shows  the  entire  6  seconds  of  the  benchmark. 
Recall  the  peaks  located  in  the  latter  half  of  the  benchmark  signal  are  a  set  of  1 1 
broadband  transients.  Notice  the  cluster  of  pulses  just  before  the  fourth  second  in  the 
graph  of  the  highest  sampling  rate.  Here,  the  three  spikes  are  easily  recognizable.  Now 
dropping  down  through  the  other  three  graphs  where  the  sampling  frequency  is  reduced, 
notice  that  the  three  spikes  are  slowly  combining  to  form  a  single  transient.  If  the  final 
stages  of  method  1  or  2  were  the  only  channels  to  pick  up  this  set  of  transients,  there 
would  be  no  way  to  differentiate  between  them. 


123 


1 

0.5 


i - 1 - 

—  f  =  1 1025  samples/sec 

— 1 - 

— r 

n 

. »i  ,1  1  1.^. 

L 

l 

2  3  4 

TIME  (seconds) 


Figure  VI-54  Reduction  of  Signal  Resolution  Due  to  Downsampling 
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VII.  CONCLUSIONS 


This  thesis  introduced  new  transient  detection  techniques  designed  to  provide 
sonar  operators  an  enhanced  capability  to  detect  enemy  acoustic  signals.  The  methods 
presented  in  this  thesis  are  based  on  the  Autoregressive-Moving  Average  and  Discrete 
Wavelet  Transform.  The  difference  is  the  order  in  which  these  two  operations  are 
performed.  Method  1  first  whitens  the  input  signal  using  an  ARMA  model  before 
denoising  the  whitened  signal  with  a  wavelet  decomposition  operation.  Method  2  first 
applies  wavelet  decomposition  before  applying  a  bank  of  whitening  filters  to  each 
wavelet  channel. 


A.  USW  LITTORAL  WARFARE  APPLICATIONS 

The  experimental  results  reveal  a  significant  increase  in  signal  detection  using 
either  proposed  method  when  the  transient  is  narrowband.  This  ability  to  localize 
frequencies  might  be  useful  for  applications  where  the  frequency  of  a  potential  transient 
is  known  a  priori,  such  as  for  example  in  active  sonar  operations  where  the  return  signal 
is  known  to  be  within  a  Doppler  shift  of  the  transmitted  frequency.  Wavelet  denoising 
has  the  ability  to  narrow  its  filtering  process  such  that  a  specific  channel  can  be  observed, 
rejecting  all  other  frequency  content.  Therefore,  this  process  has  the  potential  to  increase 
the  effective  range  of  the  sonar  without  having  to  increase  the  signal  power. 

Another  application  to  narrowband  signal  detection  might  be  for  passive  sonar, 
operations  where  one  attempts  to  uncover  a  known  frequency  within  noise  of  similar 
frequency  content.  In  low  SNR,  we  have  shown  that  the  spectrogram  does  not  provide 
useful  capabilities  in  this  regime  and  so  we  rely  on  time-domain  methods.  Examples  of 
such  signals  might  be  specific  tonals  or  blade  rates  produced  by  enemy  platforms  hidden 
within  the  noise  of  background  shipping. 

In  spite  of  the  fact  that  a  computer  can  successfully  detect  a  transient  signal  over 
some  given  threshold,  the  sonar  operator  is  still  a  vital  member  in  the  detection  phase  of 
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US W .  Indeed,  it  is  not  unusual  for  the  initial  detection  of  an  enemy  submarine  to  come 
from  a  source  other  than  an  automatic  alert  function.  Compared  to  the  other  warfare 
areas,  human  intervention  is  more  vital  in  the  detection  phase  because  the  pace  is 
generally  slower.  With  this  in  mind,  an  operator  must  have  multiple  displays  and  sensors 
at  his  disposal.  It  is  believed  that  the  methods  tested  in  this  thesis  will  not  only  increase 
detection  capabilities,  but  will  also  provide  enhanced  display  features. 

The  main  idea  behind  multiple  display  features  at  an  operator’s  disposal  is  to 
correlate  data.  Consider  a  Graphical  User  Interface  (GUI)  system  where  an  operator  is 
shown  a  wavelet  packet  framework  diagram  as  in  Figure  VI-7  where  the  boxes  represent 
the  wavelet  decomposition  channels.  The  display  feature  could  also  link  with  existing 
features  in  the  sonar  suite  to  further  enhance  transient  detection. 

Another  way  to  enhance  operator  displays  would  be  to  include  various  intensity 
levels  that  somehow  correspond  to  the  SNR.  Examples  include  full  color  displays  or 
alarms  with  a  range  of  volumes  relative  to  the  strength  of  a  transient  signal. 
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B.  SUGGESTIONS  FOR  FUTURE  RESEARCH 


As  the  wavelet  decomposition  can  theoretically  be  carried  on  to  infinite  detail,  so 
can  the  ideas  of  improving  technology.  Although  this  section  is  certainly  not  all- 
inclusive,  several  suggestions  are  provided  for  future  research  based  on  the  results  from 
Methods  1  and  2. 

1.  Method  1  Improvements 

The  DWT  denoising  process  of  Method  1  uses  the  Haar  Wavelet,  which  is  the 
most  basic  wavelet  basis  function  that  uses  only  a  first  order  polynomial  for  the  highpass 
and  lowpass  filter  coefficients.  As  one  can  see  from  the  frequency  response  of  the  filter 
coefficients  (Figure  III-9),  there  is  a  significant  overlap  between  the  two  filters,  which 
results  in  a  significant  amount  of  high  frequency  information  contained  in  the  low 
frequency  channel  and  vice  versa.  Further  research  can  be  done  to  study  other  wavelet 
basis  functions  that  have  a  sharper  frequency  response.  To  see  a  preliminary  example  of 
denoising  using  a  different  wavelet  basis  function,  refer  to  Appendix  B. 

Because  we  know  that  the  wavelet  decomposition  produces  independent  channels, 
another  area  for  follow-on  research  is  to  define  a  way  to  correlate  data  at  each  channel,  as 
described  in  Appendix  B. 

2.  Method  2  Improvements 

Recall  the  successful  improvements  the  dyadic  filtering  had  on  the  narrowband 
high  frequency  transients.  It  is  believed  that  the  same  success  will  occur  for  Method  2  if 
the  wavelet  decomposition  were  carried  out  on  the  highpass  channels  (details)  as  well. 
This  would  require  further  investigation  on  defining  the  whitening  filter  coefficients  for 
these  channels.  In  effect,  Method  2  would  then  perform  the  full  wavelet  packet 
decomposition  on  the  input  signal  and  then  run  it  through  the  modified  innovation  filter 
bank.  Like  Method  1,  the  modified  Method  2  would  result  in  15  channels  for  the  three- 
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stage  wavelet  decomposition.  Again,  this  structure  lends  itself  to  a  useful  GUI  display  in 
which  access  to  a  large  number  of  signal  data  is  relatively  simple. 

In  addition  to  investigating  the  wavelet  packet  structure  for  the  DWT-Innovation 
process,  the  same  argument  can  be  made  about  possible  improvements  using  another 
wavelet  basis  set  such  as  the  Daubechies  family.  This  process  is  a  bit  more  complicated 
than  changing  the  wavelet  set  in  Method  1  because  the  innovation  filter  bank  of  Method  2 
is  dependent  on  the  wavelet  filter  coefficients. 


3.  General  Improvements 

The  benchmark  signal  was  a  colored  noise  model  with  synthetic  signals 
representative  of  coastal  water  acoustics.  The  main  thrust  for  future  study  should  be  to 
use  actual  acoustic  data  in  this  regime.  In  this  way,  a  more  realistic  conclusion  can  be 
drawn  for  this  particular  application. 

Results  showed  that  the  denoising  process  performed  well  for  both  Methods  1  and 
2.  Perhaps  the  two  methods  could  be  combined  such  that  we  take  the  whitened  signals 
from  Method  2  and  pass  them  through  the  DWT  denoising  process  of  Method  1 . 

Looking  at  the  benchmark  signal,  note  that  the  three-second  window  performed 
well  for  the  .1  second  transient  signal.  However,  this  window  size  might  not  work  well 
for  signals  of  different  length.  If  we  wish  to  detect  a  range  of  transient  signal  durations, 
say  as  long  as  6  seconds  and  as  short  as  .01  second,  further  research  must  be  done  to 
investigate  the  effects  of  different  window  sizes  with  varying  transient  lengths. 

Finally,  the  last  suggestion  for  future  research  would  be  to  investigate  methods  for 
applying  a  neural  network  scheme  to  classify  patterns  of  signals  among  noise  in  order  to 
detect  transients  at  a  lower  SNR.  If  there  is  a  way  to  determine  features  that  positively 
identify  transient  signals,  then  for  USW  applications  it  can  be  another  valid  tool  to  make 
the  detection  process  more  robust. 
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APPENDIX  A. 


HAAR  WAVELET  PROPERTIES 
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A.  NECESSARY  TIME  DOMAIN  CONDITIONS 

This  appendix  shows  that  the  Haar  wavelet  system  meets  all  8  necessary 
conditions  for  a  valid  wavelet  system  [9].  Starting  with  the  wavelet  decomposition  tree 
in  Figure  m-6  and  the  MRA  equation,  we  wish  to  determine  the  Haar  coefficients  of  the 
PR  filter  bank,  the  scaling  function  and  the  shape  of  the  Mother  Wavelet.  Based  on  the 
set  of  eight  necessary  conditions,  we  will  show  that  the  Haar  system  of  equations  is 
indeed  a  valid  wavelet  system. 

Recall  that  all  parameters  that  define  a  wavelet  system  can  be  determined  from 
scaling  function.  If  we  can  somehow  develop  a  parameterizing  scheme  to  define  the 
scaling  coefficients  h(n) ,  we  could  determine  any  wavelet  system.  To  this  end,  we  begin 

by  choosing  an  arbitrary  filter  size  N .  The  filter  size  for  the  scaling  coefficients 
h{n)  must  be  even  and  must  satisfy  both  the  linear  constraint 


2>(n)  =  V2,  (A.l) 


N 

and  the  —  bilinear  constraints  which  leads  to 
2 


£|/z(n)|2=l.  (A.2) 


It  is  important  to  note  that  the  filter  h(n)  is  an  FIR  filter,  which  also  implies 
compact  support.  This  wavelet  property  is  highly  desirable  as  it  aids  in  the  time- 
localizing  ability  of  the  DWT  and  it  also  reduces  the  number  of  computations. 

N 

It  can  be  shown  that  there  are  —  - 1  degrees  of  freedom  (parameters)  in  choosing 

h(n )  that  will  guarantee  the  existence  of  a  valid  wavelet  system  [9].  Consider  the 
simplest  case  where  N  =  2.  This  filter  length  means  that  there  are  zero  degrees  of 
freedom  and  by  equations  (A.l)  and  (A.2),  we  get 
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fc(0)  +  Ml)  =  V2, 


(A3) 


and 


fc2(0)  +  /i2(l)  =  l. 


(A.4) 


Although  the  second  equation  is  nonlinear,  it  can  be  seen  by  inspection  that  the 
unique  solution  to  this  set  of  equations  is: 


h(n)  =  {h(0),h(l)}=lj^  J=j.  (A.5) 

These  are  the  Haar  scaling  function  coefficients,  which  were  developed  in  1910, 
separate  from  studies  in  modem  wavelet  analysis.  In  1992,  Ingrid  Daubechies  showed 
that  Haar’s  work  on  this  set  of  orthonormal  basis  functions  are  in  fact  the  simplest  set  in  a 
larger  family  of  wavelet  systems.  From  equation  (3.19),  we  can  determine  the  wavelet 
coefficients  as 

A,  (n)  =  {ft, (0), ft, (I)}=  j-i.  ^=1.  (A.6) 

With  just  the  coefficients  to  the  Haar  scaling  and  Wavelet  functions,  it  can  be 
shown  that  this  wavelet  system  meets  eight  necessary  conditions  [9]. 
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1. 


Theorem  1: 


If  <p(t)e  1}  is  a  solution  to  the  MRA  equation  (3.14),  and  if  J  cp(f)dt  &  0 ,  then 

5>(n)  =  V2.  (A.7) 

n 

For  the  Haar  Wavelet  the  above  equation  leads  to: 

h(0)+hQ)=-L+-L=j2. 

V2  v2 

2.  Theorem  2: 

If  <p(t)  is  an  I)  solution  to  the  MRA  equation  (3.14)  with  J  <p(t)dt  *  1,  then 

5>('-o=2>©=i;  (a.8) 

t  i 

with  $>(n  +  Ink)  ^  0  for  some  k,  then 

3>(2n)  =  5>(2n  +  l), 

n  n 

where  (3.30)  may  have  to  be  a  distributional  sum 
satisfied,  then  (3.30)  is  true. 

For  the  Haar  Wavelet  the  above  equation  leads  to: 

h(0)  +  h(2)  =  h(l)  +  h(3) 


(A.9) 

.  Conversely,  if  (3.31)  is 
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3. 


Theorem  3: 


-L+o=-j=+0- 

S  V2 


If  <p(t)is  an  L2  n  L1  solution  to  the  MRA  equation  (3.14)  and  if  integer  translates 
of  (p(t )  are  orthogonal  as  defined  by 

\<P<tW-W  =  ES(k)  =  ^  (A-10) 

then 

lHnMn-2^-^  (AH) 

For  the  Haar  Wavelet  the  above  equation  leads  to: 

h(0)h(-2k)  +  h(l)h(l-2k) 
h(0)h(0)  +  h(l)h(l)  =  -^ + -^  =  1  for  k  =  0 
h(0)-0  +  h(l) -0  =  0  forfc<0  and  k>0. 


Corollary  1: 

Under  the  assumptions  of  Theorem  3,  the  norm  of  h(n)  is  automatically  unity, 
i.e., 


5>(«)f=l.  (A-12) 
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For  the  Haar  Wavelet,  corollary  1  leads  to: 


/*2(0)  +  /?2(l)  =  — +  — =  1. 

2  2 


Corollary  2: 

Under  the  assumptions  of  theorem  3, 


2  h{2n)  =  Yj  H2n  + 1)  =  ~^= . 

n  n  V2 

For  the  Haar  Wavelet,  corollary  2  leads  to: 

h(0)  +  h(2)  =  h(l)  +  h(3)  =  ~^= . 

V2 


4.  Theorem  4: 

If  (pit) has  compact  support  on  0<t<N-\  and  if  (p{t-k)are  linearly 
independent,  then  h(n)  also  has  compact  support  over  0  <  n  <  N  —  1 : 

h(n)  =  0for  n  <0  and  n>  N  -1.  (A.  13) 

For  the  Haar  Wavelet  the  above  equation  leads  to: 

Haar  Wavelet  has  compact  support  as  h(n)  =  0 

for  n  <  0  and  n  >  1. 
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B.  NECESSARY  FREQUENCY  DOMAIN  CONDITIONS 


The  last  four  theorems  are  necessary  conditions  in  the  frequency  domain.  To 
begin  the  analysis  in  the  frequency  domain,  we  look  at  the  frequency  response  of  the 
scaling  coefficients  h{n) .  The  Haar  Scaling  filter  coefficients  in  the  z  -  domain  form  the 
equation 


(A.  14) 


Substituting  eJw  for  z  we  get  the  frequency  response  of  the  FIR  filter 


H0(co)=-j=( l  +  e~JW). 


(A.15) 


Likewise,  the  frequency  response  of  the  Haar  Wavelet  filter  coefficients  is 

H1(CD)  =  -j=(l-e~jw).  (A.  16) 

It  is  important  to  note  here  that  the  frequency  response  of  the  scaling  coefficients 
h(n)  (also  known  as  h0(n ))  is  that  of  a  lowpass  filter  with  its  cutoff  frequency  at  the 
Nyquist  frequency  rate  ( H  (7t)  =  0 ).  It  can  also  be  shown  that  the  frequency  response  of 
the  wavelet  coefficients  h ,  (n)  is  that  of  a  highpass  filter  with  its  maximum  at 
H(n)  =  V2  and  its  zero  magnitude  at  H( 0) .  This  is  the  property  that  allows  the  input 
signal  to  be  filtered  into  a  “constant-Q”  filter  bank  resulting  in  equal  highpass  and 
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lowpass  portions  preceding  every  downsample  as  described  in  Figure  IE-6.  Figure  IE-8 
shows  the  frequency  response  of  the  PR  filter  bank  used  for  the  Haar  Wavelet  system. 


1.  Theorem  5: 

If  <p(t)is  a  L1  solution  of  the  MRA  equation  (3.14),  then  the  following  equivalent 
conditions  must  be  true: 


5>(n)  =  tf(0)  =  V2.  (A.  17) 


For  the  Haar  Wavelet  the  above  equation  leads  to: 


H  (0)  =  (l  + 1)  =  ,  therefore 

V2 


h(0)  +  h(l)  =  H(0)  =  V2. 


2.  Theorem  6: 

For  h(n)e  1} ,  then 

'Yjhilri)  =  ^ h(2n  +  l)ifand only  if  H(n )  =  0.  (A.18) 


For  the  Haar  Wavelet  the  above  equation  leads  to: 


H  (n)  =  ~^=  (l  - 1)  =  0 ,  therefore 
\2 


h(0)  +  h(2)  =  h(l)  +  h(3)  = 
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3. 


Theorem  7: 


If  (pit)  is  a  solution  to  the  MRA  equation  (3.14)  in  L2  n  L1  and  is  a 

solution  of 


*m=T2H I 


'~W*' 


Of 

1 


(A.  19) 


and  after  iteration  becomes 


f  1  rr 

^Y| 

Wj 

mo). 


(A.20) 


then 


J  (p(i)(p{t  -  k)dt  =  S(k)  if  and  only  if  +  2  7d]2  =  1 .  (3.21) 


This  theorem  is  the  frequency  domain  equivalent  to  the  time  domain  definition  of 
the  orthogonality  of  the  scaling  function  [8].  Equation  (3.42)  shows  that  for  a  function 
that  is  narrow  in  time  is  spread  out  in  frequency. 

For  the  Haar  Wavelet  the  above  equation  leads  to: 

<£>(ty)  is  the  Fourier  Transform  of  the  scaling  function  (p(t) .  Finding  <&(co)  is  an 
iterative  process  and  equation  (3.40)  is  often  called  the  “Cascade  Theorem”.  By 
assuming  an  initial  <t>(0),  the  use  of  a  computer  to  iterate  the  process  can  show  that 
d>(<y)  converges,  its  inverse  Fourier  Transform  is  the  scaling  function  and  equation 
(3.42)  is  satisfied.  As  stated  before,  once  the  scaling  function  is  determined,  then  by 
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equation  (3.18)  one  can  determine  the  shape  of  the  Mother  Wavelet.  The  shape  of  the 
Haar  Scaling  and  corresponding  Wavelet  function  is  shown  in  Figure  TTT-7 


4.  Theorem  8: 

For  any  h{n)  e  L1 , 


y£/h(n)h(n-2k)  =  S(k)  ifandonlyif\H{cof'  +\H(o)  +  n^  =2.  (3.47) 


For  the  Haar  Wavelet  the  above  equation  leads  to: 

Starting  with  H {(o)  =  —j=  (l  +  e~JC0 ) ,  we  take  the  absolute  value  as  multiplying  by 
its  complex  conjugate  giving. 


(11 


-  +  —=re 


-JO) 


S  S  U/2  V2 


1.1  jo,')  1.1  io,  .  1  .  1 


■  +  — 7=e 


=  -  +  -ejm  +-e-ja>  +-. 
2  2  2  2 


Knowing  that  the  complex  exponential  form  of  the  cosine  function  is 


cos(x)  =  ±{eix+e-ix). 


we  substitute  this  into  the  above  expression  giving 


|#(£y)|2  =l  +  cos(ty). 


Now  taking  H{a>  +  7t)  =  -^(l  +  e  J<w+n) )  and  multiplying  by  its  complex 


conjugate  gives, 
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(J_ 

[V2 


+ 


1 

V2 


e~j«o+x) 


J 


(  ti  A  i  1  i 

_L  +  _l=eKco+’I)  =  -  +  -ejtoejn  +-e~JweJX 

U/2  &  2  2  2 


Again  we  substitute  the  cosine  function  in  for  its  complex  exponential  form 

except  now  it  is  shifted  by  180  degrees  because  of  Euler’s  Identity  e  *  =-l.  The 
resulting  expression  yields 


|#(<y+;r)|2  =l-cos(<y). 


Now  we  add  both  expressions  together  resulting  in  the  final  expression 
|#(<y)|2  +  \h(co  +  =  1  +  cos(ty)  + 1  -  cos (co)  =  2 . 
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APPENDIX  B.  FURTHER  RESEARCH 


The  denoising  procedure  used  in  Method  1  uses  the  Haar  Wavelet  system,  which 
is  the  most  basic  wavelet  set.  Looking  at  the  wavelet  decomposition  as  an  inner  product 
between  the  basis  set  and  the  signal,  we  see  that  the  Haar  wavelet  (Figure  III-8)  works 
well  with  functions  that  are  piecewise  continuous,  but  not  smooth  as  in  the  case  of  most 
signals.  Therefore  it  would  be  better  to  find  a  basis  set  where  the  Mother  Wavelet  has  a 
shape  similar  to  that  of  the  signal,  such  as  a  wavelet  from  the  Daubechies  family  [15] 
shown  in  Figure  B-l.  By  denoising  using  a  smoother,  compact  wavelet  we  still  get  the 
highpass  and  lowpass  effect,  but  now  the  wavelet  coefficients  (output  of  the  filter)  might 
be  larger  due  to  better  correlation  between  the  signal  and  wavelet.  Consider  the  Haar 
Wavelet  in  the  frequency  domain  shown  in  Figure  1H-9.  It  has  been  shown  that  the  Haar 
filter  coefficients  are  only  first-order  polynomials.  These  basic  filter  coefficients  imply 
that  the  cutoff  frequencies  of  the  highpass  and  lowpass  channels  overlap  significantly. 
Therefore,  the  DWT  denoising  process  might  perform  better  with  a  filter  pair  whose 
frequency  response  is  sharper  at  the  cutoff  frequency  like  Daubechies  10  (dblO)  as  shown 
in  figure  B-2.  Indeed,  upon  one  simple  test  using  the  same  5  kHz  transient  signal,  the 
dblO  wavelet  system  was  able  to  reduce  the  noise  threshold  lower  than  the  Haar  wavelet, 
as  shown  in  Figures  B-3  through  B-6. 
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Figure  B-l  Daubechies  10  (dblO)  Wavelet 


Figure  B-2  Frequency  Response  of  dblO  Wavelet  System  (fs=11.025  kHz) 


Figure  B-3  Wavelet  Packet  Framework  (AAD3  and  DAD3 ) 


LOWPASS  FILTER  (APPROXIMATION) 
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Figure  B-4  Denoising  5  kHz  Signal  Using  Haar  Wavelet  System  ( AAD3  and  DAD3 ) 


Figure  B-5  Wavelet  Packet  Framework  ( AAD3  and  DAD, ) 


Figure  B-6  Denoising  5  kHz  Signal  Using  dblO  Wavelet  System  (AAD3  and  DAD3 ) 
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Another  modification  to  Method  1  that  might  be  considered  for  further  research  is 
to  define  a  Figure  of  Merit  (FOM),  which  would  allow  for  comparison  of  the  various 
wavelet  channels.  To  this  end,  we  investigate  the  statistics  of  white  noise.  By  definition, 
a  zero  mean  white  noise  process  x(t )  is  defined  by  its  first  moment  (or  mean)  described 

as 


£{x(r)}=0  (B.l) 

and  second  moment  (or  variance)  described  as 

E{x(t)x(t  -  t)}  =  <tw2S(t)  .  (B.2) 

In  general,  we  wish  to  come  up  with  an  operation  that  can  describe  the  correlation 
between  N  white  noise  processes,  where  N  is  the  number  of  channels  in  the  wavelet 
decomposition  and  the  resulting  correlation  matrix  Rxx  is  of  dimension  N  x  N .  In 

addition,  this  operation  must  provide  visual  clarity  in  distinguishing  a  transient  signal 
from  background  noise.  With  this  in  mind,  a  figure  of  merit  is  developed  as 


FOM  =  1- 


mingig(£^) 
V  rna xeig(R„) 


(B.3) 


The  ratio  of  the  minimum  eigenvalue  to  the  maximum  eigenvalue  will  always  be 
a  number  between  zero  and  unity.  If  the  minimum  eigenvalue  is  equal  to  the  maximum, 
indicating  a  set  of  completely  uncorrelated  signals  with  normalized  power,  then  the  ratio 
would  equal  1  and  the  whole  FOM  would  equal  zero.  On  the  other  end,  if  there  is  some 
correlation  between  any  two  signals  (i.e.  a  transient  detected  in  at  least  one  channel  of 
the  multiresolution  filter  bank),  the  ratio  in  the  expression  (B.3)  will  be  less  than  unity 
and  the  entire  expression  will  be  greater  than  zero.  This  non-zero  value  will  show  up  as  a 
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“spike”  when  the  FOM  is  plotted  as  a  function  of  time.  The  square  root  simply  smoothes 
out  the  FOM  so  that  an  indication  of  a  transient  signal  is  easier  to  see  on  a  visual  display. 

There  is  one  final  note  about  this  method.  Because  of  the  downsampling 
operation  after  the  filter  bank,  the  length  of  data  will  decrease  by  a  factor  of  two  at  each 
stage.  In  order  to  make  every  vector  the  same  length,  the  signal  is  upsampled  and  run 
through  an  averaging  filter.  The  effect  of  this  process  “stretches”  the  signal  to  the  desired 
number  of  samples,  and  makes  every  other  sample  the  average  of  its  two  neighbors. 
Performing  this  operation  does  not  affect  the  statistics  of  the  signal,  nor  does  it  wash  out 
any  transient  information.  The  other  option  to  perform  the  correlation  operation  is  to 
downsample  every  detail  to  the  size  of  the  finest  approximation.  This  too  would  be  a 
valid  process,  but  with  downsampling  comes  the  loss  of  information.  The  other  reason  to 
upsample  and  average  is  that  the  correlation  function  determined  by  a  computer  is  an 
averaging  operation.  The  upsample-and-average  method  provides  more  data  points  to 
average  and  therefore  a  better  correlation  representation.  After  calculating  the  details  and 
approximations,  the  data  is  run  through  an  upsample/averaging  process  to  achieve  the 
proper  data  length  (indicated  as  primed  details  and  approximations  in  Figure  B-7). 

The  re-sizing  of  the  detail  and  approximation  at  the  finest  scale  are  then  combined 
together  as  a  set  of  column  vectors  as  shown  in  Figure  B-8  and  in  the  equation 


w  ■ 


t  t 

dl  d2 

i  i 


T 

d . 


(B.4) 


146 


147 


For  this  procedure,  three  stages  of  multiresolution  are  again  used.  Finally,  the 
correlation  function  is  computed  as 


R 


WW 


w, 


(B.5) 


and  the  figure  of  merit  is  determined  as 


FOM  (f)  =  l- 


mmeig 


ij  max  eig(R„„(t)) 


(B.6) 


A  plot  of  the  Figure  of  merit  reveals  values  going  from  zero  to  unity  as  correlation 
between  any  two  inputs  increases.  A  “spike”  on  the  FOM  plot  reveals  a  transient  signal. 
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APPENDIX  C.  MATLAB  PROGRAMS 
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A.  CODE  COMMON  TO  BOTH  METHODS 


1.  ARMA  model 

function  [b, a] =cristi_ar (y) ; 

% 

%  [b,a] =cristi_ar (y) ; 

a 

o 

%  y  =  data  (a  COLUMN  vector) 

ry=cristi__corr  (y,y)  ; 

%  Levinson  recursion 

a=l ;  b=l ;  sig2=ry(l);  sig2p=sig2; 

p-0 ;  fact=0; 

%  while  (p<10) ; 

while  (p<10)&(fact<0.98) ; 

r=ry (2  :p+2)  ; 

D=r ' *f lipud(a) ;  Dprime=r ' *f lipud (b) ; 

g=D/sig2p;  gprime=Dprime/sig2 ; 

ap= [a; 0] -g* [ 0 ; f lipud (b) ] ;  b= [b; 0] -gprime* [0 ; f lipud (a) ] ; 
a=ap; 

f act=l-g*gprime ; 
sig2=fact*sig2 ; 
sig2p=fact*sig2p; 

P=P+1; 

end 

b=sqrt (sig2p) ; 


2.  Autocorrelation  function 

function  Rxy=cristi_corr (x,y) 

% Function  will  calculate  the  cross-correlation  function. 
%It  is  faster  than  MATLAB  function  XCORR 
%  Rxy=cristi_corr (x, y) 

% 

m=rnax(  [ length  (x)  ,  length  (y)  ] )  ; 
fx=f  ft  (x,  2*m)  ; 
fy=f  ft  (y,2*m)  ; 

Sxy=fx. *conj (fy) ; 

Rxy=real (if ft (Sxy) ) /m; 

%end  %  cristi__corr 
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3. 


Benchmark  Test 


%Name:  John  D.  Stevens 

%Date :  22  Jul  1999 

%Filename:  meth_0102.m 

%description :  this. program  uses  Methods  1  and  2  to  detect 

various 

%transient  signals  in  a  colored  noise  environment. 

% 

clear 

%  [s,  fs,bps]  =wavread(  'H:  \matlab2\thesis\data\wav_data\record' )  ; 

% school 

% [ S , f  s , bps ] =wavread ( ' C :  \My  • 

Doctiments\dad\m_f iles\thesis\data\wav_data\record' ) ;  %home 
%load  'H:\matlab2\thesis\data\mat_data\fan.mat' 
load  'H: \matlab2\thesis\data\mat_data\fan2 .mat ' 

%load  'C:\My  Documents\dad\m_files\thesis\data\wav_data\fan.mat' 
%load  'C: \My  Documents \dad\m_f iles\thesis\data\wav_data\fan2 .mat ' 
fs=11025;  %telephone  quality 

^initialize  the  confidence  vector: 

%narrowband  signals: 

nbl= [50; .5; .5] ;nblh= [50; . 1 ;  .5]  ; 

nb2  =[200; .1; .1] ;nb2h=[200;  .2;  .2] ; 

nb3= [600; .1; .05] ;nb3h=[600;  .1;  .08]  ; 

nb4= [900; .1; .03] ;nb4h=[900; .1;  .05]  ; 

nb5= [2000 ; .1; .01] ;nb5h=[4000;  .1; .1]  ; 

nb6= [5000 ; .1; .01] ;nb6h=[5000; .1; .003] ; 

NB= [nbl , nb2 , nb3 , nb4 , nb5 , nb6 ] ; 

NBH=[nblh,nb2h,nb3h,nb4h/nb5h,nb6h] ; 

BBC=[ .1; .8] ; 

BBW=[ .1;  .01]  ; 
trials=length(NB); 
max_dtr=max (NB ' ) ; 
max_dtr=f loor ( fs*max_dtr (2 ) )  ; 

fig=l; 

%TRSOUND= [ ] ; 

%CSOUND= [ ] ; 

%WSOUND= [ ] ; 

%for  c=l:trials; 
c=6; 

%ftr=NB(l,c) ; 
ftr=NBH(l,c) ; 

%dtr=NB(2, c) ; 
dtr=NBH (2 , c) ; 

%dtr=BBC (1) ;  %duration  of  colored  noise  transient 
%dtr=BBW(l) ;  %duration  of  white  noise  transient 

%atr=NB (3 , c) ; 
atr=NBH (3 , c) ; 
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%atr=BBC (2) ; 

%atr=BBW (2 ) ; 

%create  transient  signal: 

%A)  narrowband  sinusoidal  transient 

NS_T=floor (fs*dtr) ;  %number  of  samples  in  the  transient 

tot__time=NS_T/f  s ; 

t=0:  (tot_time/  (NS_T-1)  )  :tot_time; 

y=cos (2*pi*f tr*t ) ; 

ham_win=hamining  ( length  ( t )  )  ; 

y=y .  *ham_win '  ; 

%B)  broadband  transient  signal: 

%NS_T= floor ( fs*dtr) ;  %number  of  samples  in  the  transient 
%y=randn(l,NS_T) ; 

%ham_win=hamming (NS_T) ; 

%y=y. *ham_jwin' ; 

%C)  chirp  signal: 

%NS_T=floor (fs*dtr) ;  %number  of  samples  in  the  transient 
%tot__time=NS_T/f  s; 

%t=0:  (tot_time/  (NS_T-1)  )  :tot_time; 

%tp= [0  dtr/4  dtr/2  3*dtr/4  dtr] ;  %  time  breakpoints 
%fp=[0  200  100  150  300] ;  %  instantaneous  frequency 

breakpoints 

%p=polyf it ( tp, fp/ 4) ;  %  fit  4th  order  polynomial  over  time 

%y=chirp ( t ,  p) ; 

%D)  colored  noise  transient  signal: 

%NS_T=f loor (f s*dtr) ;  %number  of  samples  in  the  transient 
%y=fan2 (20000 : 20000+NS_T-l) ' ; 

%ham_ win=hamming  (NS_T)  ; 

%y=y .  *ham_win '  ; 

%E)  colored  noise  transient  signal:  ELVIS 
%load  ' C : \My 

Documents\dad\m_f iles\ thesis \data\wav_data\elvis_a .mat # 

%load  'C:\My 

Document  s\dad\m__files\thesis\  data  \wav_data\elvis__b.  mat ' 

%load  '  H: \matlab2\ thesis \data\ma t„data\elvis_a. mat ' 

%load  ' H : \matlab2 \ thesis \data\ma t_data\elvis_b .mat 9 
%NS_T=length(elvis_a) ; 

%y=elvis_a ' ; 

% hamjwin = hamming  ( N S_T )  ; 

%y=y  *  *ham_win  '  ; 

%atr= . 3  ? 

yf=fan2; 

yf  (10001:  ( 10000 +NS_T) ) =yf (10001:  (10000+NS_T) )+atr*y'; 

%F)  use  the  data  fan  (includes  the  tapping)  as  the  input 

signal 

%yf =f an; 

%method  1:  pre-whitening 
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% [zO,  zll_a, zlh_a, z21_a, z2h_a, z31_a, z3h_a] =method_01_f ilt (yf ) ; 

% [zO, zll_a, zlh_a, z21_a,  z2h_a, z31_a, z3h_a, SKR_1] =method_01_f ilt 
v(yf ) ; 

%  [zO, zll_a, zlh_a, z21_a, z2h_a, z31_a, z3h_a] =tree_filt (yf ) ; 

[ STG1 , STG2 , STG3 , STG4 ] =wavepackz (yf ) ; 

%  SNR_1 =my_SNRl ( STG1 , STG2 , STG3 , STG4 ) ; 

%  SNR1 =my_SNRl 4p 1 o t s (STG1, STG2 , STG3 , STG4) ; 

%method  2:  post -whitening  (Cristi's  algl) 

[zll_b, zlh_b, z21_b, z2h_b, z31_b, z3h_b] =method_02_filt (yf )  ; 

% [zll_b, zlh_b, z21_b, z2h_b, z31_b, z3h_b, SNR_2] =method_02_f iltv(y 
f)  ; 

%SNR2=my_SNR24plots (yf ) ; 

%method  3:  SVD  model  (Cristi's  alg2) 

% [zl_c, z2_c, z3_c, z4_c, z5_c( z6_c] =method_03_f ilt (yf ) ; 

%y= [y, zeros (l,max_dtr-length(y) ) ] ; 

%TRSOUND=[TRSOUND;y] ; 

%CSOUND=[CSOUND;yf ' ] ; 

%WSOUND= [WSOUND; zO ' ] ; 

%plot  results: 
s3_a=3 ; i3_a=3+s3_a; 
s3_b=4 ; i3_b=3+s3_b; 
s4_a=5; i4_a=7+s4_a; 
s4_b=6 ; i4_b=7+s4_b; 

fig=plot_meth0102 (STG1, STG2 (1, : ) ,STG2(2, : ) , STG3 (s3_a, : ) , STG3 (s 
3_b, :) _ 

STG4 (s4_a, : ) , STG4 (s4_b, : ) ,NS_T, fig) ; 

%f ig=plot_meth0102 (zO, zll_a, zlh_a, z21_a, z2h_a, z31_a, z3h_a,NS_T 
/ fig) ; 

f ig=plot_meth0102 (STG1, zll_b, zlh_b, z21_b, z2h_b, z31_b, z3h_b,NS_ 
T, fig) ; 

%f ig=plot_meth03 (zO , zl_c, z2_c, z3_c , z4_c, z5_c, z6_c,NS_T, fig)  ; 
figure (fig) ; 

specgram(yf , 512 , fs, 512 , 400) ; 
fig=fig+l; 

%end 
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4.  Plotting  Routines 


%Name:  John  D.  Stevens 

%Date :  26  Feb  2000 

%Filename:  plot _jneth0102  .m 

%description:  this  function  will  run  a  signal  through  the  first 

algorithm 

%that  was  developed  by  Dr.  Roberto  Cristi  and  plot  it 
% 

% 

b=plot_meth0102 ( zf 0 , zf Ot , zf 1 , zf It , zf 2 , zf 2t , zfar , NS_T, fig_num) 

% 

function 

b=plot_meth0102 (zfar, zf 0 , zf Ot , zf 1, zf It , zf 2 # zf2t , NS_T, fig_num) ; 
TTL= [  "  ]  ; 
slp=l ; 

%fig__num=l  ; 

% indexes  due  to  downsampling: 

indl=10000 ; 

ind2=f loor ( indl/2 ) ; 

ind4  =  floor (indl/4)  ; 

ind8  =  floor (indl/8)  ; 

Ns_tl=NS__T; 

Ns_t2=f loor (NS_T/2 ) ; 

Ns_t4=f loor (NS_T/4) ; 

Ns_t8  =  f  loor  (NS__T/8 )  ; 

%create  the  spikes  surrounding  the  signal 
sfar= zeros (size (zfar) )  ; 
sfar (indl) =max(zfar) ; 

sfar (indl+floor (Ns_tl*slp) ) =max(zfar) ; 

sf0=zeros (size(zfO) ) ; 
sfO (ind2) =max(zf0) ; 

sfO (ind2+floor (Ns_t2*slp) )=max(zf0) ; 

sf 0t=zeros (size (zf Ot) ) ; 
sfOt (ind2) =max(zf Ot) ; 

sf  Ot  (ind2+f  loor  (Ns__t2*slp)  )  =max(zf0t)  ; 

sfl=zeros (size (zfl) ) ; 
sfl (ind4) =max(zfl) ; 

sfl (ind4+floor (Ns_t4*slp) ) =max(zfl) ; 

sflt=zeros (size(zflt)); 
sflt (ind4) =max (zf It) ; 

sf 1 t { ind4+ floor (Ns_t4*slp) ) =max(zflt) ; 

sf2=zeros (size (zf 2) ) ; 
sf2 (ind8) =max(zf2) ; 
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sf2  (ind8+f loor (Ns_t8*slp)  )  =xnax(zf2)  ; 

sf2t=zeros (size(zf2t) ) ; 
sf2t (ind8) =max(zf2t) ; 

sf2t (ind8+f loor (Ns_t8*slp) )=max(zf2t) ; 

%create  the  time  vectors: 

Nar=length{zfar) ; 

NO=length(zfO) ; 

Nl=length(zfl) ; 

N2=length{zf2) ; 

fs=11025; 
f_t ime=Nar / f s ; 

t_ar=0 : f_time/ (Nar-1) :f_time; 
t_0=0 : f_time/ (NO-1) :f_time; 
t__l=0  :  f_time/  (Nl-1)  :f_time; 
t_2  =  0:f_time/  (N2-1)  :f_time; 

%plot  the  figures 
f igure (f ig_num) 

plot  (t_ar,  zfar,  'b'  ,-t_ar,  sfar,  'r :  ' )  ; 
title ( [TTL, 'AR  WHITENING  FILTER' ] ) ; 
ylabel ( ' MAGNITUDE ' ) ; 

axis(  [0,  f_time/min(zfar)  /iaaxfzfar)  ]  )  ; 
xlabel { ' TIME  ( seconds )  ' )  ; 

%ylabel ( 'MAGNITUDE' ) ; 

%legend( [ 'SNR  =  ' , num2str ( SNR (1) ) , '  dB']); 
f i gur e ( f i g_num+ 1 ) 

subplot  (2,1,1)  ,.plot  ( t— 0 ,  zf  0,  'b'  ,  t_0 ,  sf  0 ,  '  r :  ' )  ; 

%title( [TTL, 'LOWPASS  FILTER  STAGE 
1' ] ) ,axis ( [0, f_time,min(zfO) ,max(zf0) ] ) ; 
title ( [TTL, 'LOWPASS  FILTER 
(APPROXIMATION) ' ] ) ,axis ( [ 0 , f_time , 0 , 1 ] ) ; 

%ylabel ( 'MAGNITUDE' ) ; 

%Iegend( [ 'SNR  =  ' , num2str (SNR(2 ) ) , '  dB']); 
subplot (2,1,2) ,plot ( t_0 , zf Ot , 'b' , t_0 , sf Ot , ' r : ' ) ; 

%title( 'HIGHPASS  FILTER  STAGE 

1' ) ,axis ( [0, f_time,min(zf Ot) ,max(zf0t) ] ) ; 

title ( [TTL, 'HIGHPASS  FILTER  (DETAIL) ']), axis ([ 0 , f_time , 0 , 1 ] ) 
xlabel ('TIME  (seconds)'); 

%ylabel { ' MAGNITUDE ' ) ; 

%legend{ [ 'SNR  =  ' , num2 str ( SNR ( 3 ) ) , '  dB ' ] ) ; 
figure ( f ig_num+2 ) 

subplot (2,1,1) ,plot (t_l, zfl, 'b' , t_l, s  f 1 , ' r : ' ) ; 

%title( [TTL, 'LOWPASS  FILTER  STAGE 
2 ' ] ) ,axis ( [0, f_time,min (zf 1) ,max(zfl) ] ) ; 
title ( [TTL, 'LOWPASS  FILTER 
(APPROXIMATION) '] ) ,axis( [ 0 , f_time, 0 , 1 ] ) ; 

%ylabel ( 'MAGNITUDE ' ) ; 

%legend( [ 'SNR  =  ' , num2str (SNR(4) ) , '  dB']); 
subplot  (2,1,2)  ,plot  ( t_JL,  zf  It ,  'b'  ,  t_l ,  sf  It ,  '  r :  ' )  ; 
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%title( 'HIGHPASS  FILTER  STAGE 
2 ' ) , axis ( [ 0 , f_time,min (zf It ) ,max (zf It) ] ) ; 

title ( [TTL, 'HIGHPASS  FILTER  (DETAIL) ']), axis ([ 0 , f_time, 0 , 1 ]) ; 
xlabel ( 'TIME  (seconds) ' ) ; 

%ylabel ( 'MAGNITUDE ' ) ; 

%legend( [ 'SNR  =  ' ,num2str (SNR(5) ) , '  dB']); 
figure ( f ig_num+3 ) 

subplot (2 , 1,1) , plot (t_2 , zf2 ( 'b' ,t_2,sf2, 'r: ') ; 

%title( [TTL, 'LOWPASS  FILTER  STAGE 
3'  ]  )  , axis(  [0,  f_time,min(zf2)  ,max(zf2)  ]  )  ; 
title ( [TTL, 'LOWPASS  FILTER 
(APPROXIMATION) ' ] ) , axis ( [0 , f_time, 0,1]); 

%ylabel ( 'MAGNITUDE '); 

%legend( [ 'SNR  =  ' ,num2str(SNR(6) ) , '  dB']); 
subplot (2,1,2) ,plot (t_2, zf2t, 'b' ,t_2,sf2t, 'r: ') ; 

%title( 'HIGHPASS  FILTER  STAGE 
3 ' ) ,axis ( [0 , f_time, min ( zf2t) ,max(zf2t) ] ) ; 

title ( [TTL, 'HIGHPASS  FILTER  (DETAIL) ']), axis ([ 0 , f_time, 0 , 1] ) ; 
xlabel ( 'TIME  (seconds) ' ) ; 

%ylabel ( 'MAGNITUDE ' ) ; 

%legend( [ 'SNR  =  ' ,num2str(SNR(7) ) , '  dB']); 
b=fig_num+4; 


%Name:  John  D.  Stevens 

%Date :  22  Jul  1999 

%Filename:  SNR_plots3.m 

%description:  plots  the  results  from  my_SNRl  and  my_SNR2 


clear 

%load  up  results: 

load  'H:\matlab2\thesis\data\mat_data\SNR_ARb.mat' 
load  'H: \matlab2\thesis\data\mat_data\SNR_Mlb.mat' 
load  ' H : \matlab2 \ thesis \data\mat_data\ SNR_M2b . mat ' 
%load  'C:\My 

Documents\dad\m_files\thesis\data\wav_data\snr_ar.mat' 

%load  'C:\My 

Documents\dad\m_files\thesis\data\wav_data\snr_Ml.mat' 

%load  'C:\My 

Documents \dad\m_f iles\thesis\data\wav_data\snr_M2 .mat ' 

FTR=0 =10:5000 ; 

ATR=0 : .002:1; 

thresh=3 ;  %threshold  in  dB 

ax= [min (FTR) ,max(FTR) ,min(ATR) , max (ATR) ,-20,50]  ; 


figured) 

%surf (FTR, ATR, SNR_AR) ;axis (ax) ; 
mesh ( FTR , ATR , SNR_AR , SNR_M1 ) ; axi s(ax); 

C=colormap; 

title ( 'NARROWBTUgD  DETECTION  REGION  FOR  WHITENING  FILTER'); 
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xlabel (' FREQUENCY  OF  TRANSIENT  (Hz)'); 
ylabel ( 'AMPLITUDE  OF  TRANSIENT' ) ; 
zlabel ( ' SNR  (dB)'); 

figure (2) 

%  sur  f ( FTR , ATR , SNR_M1 ) ; axi s ( ax ) ; 
mesh ( FTR , ATR , SNR_M1 , SNR_M1 ) ; axi s ( ax ) ; 

title ( 'NARROWBAND  DETECTION  USING  WAVELET  DENOISING' ) ; 
xlabel (' FREQUENCY  OF  TRANSIENT  (Hz)'); 
ylabel ( 'AMPLITUDE  OF  TRANSIENT ' ) ; 
zlabel ('SNR  (dB) ' ) ; 

figure (3) 

%surf  ( FTR,  A.TR,  SNR_M2 )  ; axis  (ax)  ; 
mesh ( FTR , ATR , SNR_M2 , SNR_Ml ) ; axi s ( ax ) ; 

title ( 'NARROWBAND  DETECTION  USING  MULTIRESOLUTION  INNOVATION  - 
THRESHOLD  =  3  dB ' ) ; 

xlabel ( 'FREQUENCY  OF  TRANSIENT  (Hz) ' ) ; 
ylabel ( 'AMPLITUDE  OF  TRANSIENT' ) ; 
zlabel ('SNR  (dB) ' ) ; 

figure (4) 
x=SNR_AR; 

I=find(x<thresh)  ; 
x(I) =thresh; 

contourf (FTR, ATR,x) ;colorbar; 
colormap ( jet)  ; 

%pcolor ( FTR , ATR , x ) ; 

title ( 'NARROWBAND  DETECTION  REGION  FOR  WHITENING  FILTER  - 
THRESHOLD  =  3  dB ' ) ; 
xlabel ( ' FREQUENCY  (Hz ) ' ) ; 
ylabel ( 'AMPLITUDE ' ) ; 

figure (5) 
x=SNR_Ml ; 

I=find(x<thresh) ; 
x(I) =thresh; 

contourf (FTR, ATR, x) ;colorbar; 
colormap (jet) ; 

%pcolor (FTR, ATR, x) ; 

tit le( 'NARROWBAND  DETECTION  USING  WAVELET  DENOISING  -  THRESHOLD 
3  dB ' ) ; 

xlabel ( ' FREQUENCY  ( Hz ) ' ) ; 
ylabel ( 'AMPLITUDE ' ) ; 

figure (6) 
x=SNR_M2 ; 

I=find(x<thresh) ; 
x(I)=thresh; 

contourf (FTR, ATR, x) ;colorbar 
colormap (jet) ; 

%pcolor (FTR, ATR,x) ; 

title ( 'NARROWBAND  DETECTION  USING  MULTIRESOLUTION  INNOVATION  - 
THRESHOLD  =  3  dB '  )  ; 
xlabel ( ' FREQUENCY  ( Hz ) ' ) ; 
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ylabel ( 'AMPLITUDE' ) ; 

%DETERMINE  PERCENT  IMPROVEMENT  ON  DETECTION: 

%A)  FOR  THE  ENTIRE  SPECTRUM: 
tot_pts_AR=prod ( size ( SNR_AR) ) ; 

d.et_pts_AR=find(SNR_AR>=thresh)  ;  %find  all  points  >=  3db; 
det_pts_AR=length (det_pts_AR) ; 

pct_det_AR= (det_pts_AR/tot_pts_AR) *100 ;  %percent  detection  for 
whitening  filter 

tot_pts_Ml=prod(size (SNR_M1) ) ; 

det_pts_Ml=find(SNR_Ml>=thresh) ;  %find  all  points  >=  3db; 
det_pts_Ml=length(det_pts_Ml) ; 

pct_det_Ml= (det_pts_Ml/tot_pts_Ml) *100;  %percent  detection  for 
whitening  filter 

pct_incrMl=100* (pct_det_Ml-pct_det_AR) /pct_det_AR;  %detection 
percent  increase 

tot_pts_M2=prod(size (SNR_M2) ) ; 

det_pts_M2=f ind (SNR_M2>=thresh) ;  %find  all  points  >=  3db; 
det_pts_M2=length(det_pts_M2 ) ; 

pct_det_M2= (det_pts_M2/tot_pts_M2 ) *100 ;  %percent  detection  for 
whitening  filter 

pc t_inc rM2 =100* (pct_det_M2-pct_det_AR)  /pct_det_AR;  %detection 
percent  increase 

PCT_DET_TOT= [pct_det_AR;pct_det_Ml ; pct_det_M2  ] 

PCT_BETTER= [ 0 ; pct_incrMl ;pct_incrM2 ] 

%B)  FOR  LOW  FREQUENCIES  BELOW  2kHZ  ONLY: 
f in=find (FTR==2000) ;  %column  where  freq‘=  2000  Hz 

tot_2k_AR=prod (size ( SNR_AR ( : , 1 : fin) ) ) ; 

det_2k_AR=f ind ( SNR_AR ( : , 1 : fin) >=thresh) ;  %find  all  points  >= 
3db; 

det_2k_AR=length(det_2k_AR) ; 

pct_2k_AR= (det_2k_AR/tot_2k_AR) *100 ;  %percent  detection  for 
whitening  filter 

tot_2k_Ml=prod ( size ( SNR_M1 ( : , 1: f in) ) ) ; 

det_2k_Ml=f ind(SNR_Ml ( : , 1 : fin) >=thresh) ;  %find  all  points  >= 
3db; 

det_2k_Ml=length(det_2k_Ml) ; 

pct_2k_Ml= (det_2k_Ml/tot_2k_Ml) *100 ;  %percent  detection  for 
whitening  filter 

pc t_2 k_inc rMl =100* (pct_2k_Ml-pct_2k_AR) /pct_2k_AR;  %detection 
percent  increase 

tot_2k_M2=prod ( size ( SNR_M2 ( : ,l:fin) ) ) 

det_2k_M2=f ind (SNR_M2 ( : , 1 : fin) >=thresh) ;  %find  all  points  >= 
3db; 

det_2k_M2=length(det_2k_M2) ; 

pct_2k_M2= (det_2k_M2/tot_2k_M2) *100;  %percent  detection  for 
whitening  filter 
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pct_2k_incrM2=100* (pct_2k_M2-pct_2k_AR) /pct_2k_AR;  %detection 
percent  increase 

PCT_DET_2k= [pct_2k_AR;pct_2k_Ml ;pct_2k_M2 ] 

PCT_BETTER_2k= [ 0 ; pct_2k_incrMl ; pct_2k_incrM2 ] 


%C)  FOR  LOW  FREQUENCIES  BELOW  300  Hz  ONLY: 
fin=find(FTR==300) ;  %column  where  freq  =  300  Hz 

tot_3c_AR=prod(size (SNR_AR( : , 1 : fin) )); 

det_3c_AR=find(SNR_AR( : , 1 : f in) >=thresh) ;  %find  all  points  >= 
3db; 

det_3c_AR= length (det_3c_AR) ; 

pct_3c_AR= (det_3c_AR/tot_3c_AR) *100;  %percent  detection  for 
whitening  filter  ■> 

tot_3c_Ml=prod(size (SNR_M1 ( : , 1 : fin) ) ) ; 

det_3c_Ml=find(SNR_Ml ( : , 1 : fin) >=thresh) ;  %find  all  points  >= 
3db; 

det_3c_Ml=length(det_3c_Ml) ; 

pct_3c_Ml= (det_3c_Ml/tot_3c_Ml) *100;  %percent  detection  for 
•whitening  filter 

pc t_3 c_inc rMl = 1 0 0 * (pct_3c_Ml-pct_3c_AR) /pct_3c_AR;  %detection 
percent  increase 

tot_3c_M2=prod{size(SNR_M2 ( : , l:fin) ) ) ; 

det_3c_M2=find(SNR_M2 ( : , l:fin) >=thresh) ;  %find  all  points  >= 
3db; 

det_3c_M2=length(det_3c_M2) ; 

pct_3c_M2= (det_3c_M2/tot_3c_M2) *100;  %percent  detection  for 
whitening  filter 

pc t_3 c_incrH2  =  100* (pct_3c_M2-pct_3c_AR) /pct_3c_AR;  %detection 
percent  increase 

PCT_DET_3 c= [ pc  t_3  C_AR ; pc  t_3  C_M1 ; pc  t_3  C_M2 ] 

PCT_BETTER_3c= [ 0 ;pct_3c_incrMl ;pct_3c_incrM2 ] 
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5.  Display  Filter 


%Name:  John  D.  Stevens 

%Date :  21  Jan  2000 

%Filename :  smoother . m 

%description :  this  program  will  take  a  white  noise  process  and 

output 

%a  normalized  smoothed  version  for  display.  This  process  also 
%gives  an  approximation  to  the  noise  variance  and  the  signal 
variance, 

%which  is  useful  to  compute  the  signal  to  noise  ratio 
% 

%  zf=smo other (y) 

% 

function  zf=smoother (y) 

if  size (y, 2 ) >size (y, 1) ,  y=y' ;  end  %  make  y  into  column  vectors 

z=y. A2; 

B=1 ; 

A= [  1 ,  -  .  9  9  ]  ; 

zf=f ilter (1, [ 1 , - . 99 ] ,z) ; 

zf =zf ' ;  %make  output  in  row  vectors 

^normalize  the  results: 

mzf=max(zf ,  [] ,2) ; 

zf=zf . /mzf { : , ones (length (zf)  ,  1)  )  ; 

%  [zf ,  f_f in]  =f ilter  (B,  A,  z ,y_init )  ; 


6.  Broadband  Transient  SNR  Test 

%Name:  John  D.  Stevens 

%Date :  22  Jul  1999 

%Filename :  SNR„elvis 

%description: 


clear 

load  'H : \matlab2\thesis\data\mat_data\fan2 .mat # 
load  'H: \matlab2\ thesis \data\ma t_data\elvis_a .mat ' 

%load  ' C : \My  Documents \dad\m_f iles\thesis\data\wav_data\f an2 . mat # 
%load  'C:\My 

Document s  \ dad \m_f  i  1  e s  \ the s i s \ da t a  \ wav_„da t a  \  e  1  vi s_a .  mat ' 

f s=11025 ;  %telephone  quality 

SNR_AR= [ ] ; 

SNR__M1  =  [  ]  ; 

SNR_M2=[]; 

dtr= . 1 


160 


%NS__T=length  (elvis__a)  ; 

NS_T=floor (fs*dtr) ;  %number  of  samples  in  the  transient 
ham_win=hamming  (NS_T)  ; 

% ELVIS 
%y=elvis_a ' ; 

%y=y/max(y);  %normalize  the  elvis  signal 
%y=y .  *ham__win '  ; 

% WHITE  NOISE: 
y=randn(l,NSJT) ; 
y=y . *ham_win ' ; 

for  atr=0 : . 005 : 1 ;  %amplitude 
atr 

%add  transient  to  benchmark  signal 
yf=fan2 ; 

yf  (10001:  (10000+NS_T) )=yf (10001: (10000+NS_T) )+atr*y' ; 

%method  1:  pre -whitening 
[ STG1 , STG2 , STG3 , STG4 ] =wavepack  (yf )  ; 
snr_l=my__SNRl ( STG1 , STG2 , STG3 , STG4 ) ; 

%method  2 :  post-whitening 
snr_2  =my_SNR2  (yf )  ; 

SNR_AR= [ SNR_AR ; snr_l ( 1 )  ]  ; 

SNR_M1= [ SNR_M1 ; snr_l ( 2 ) ] ; 

SNR_M2= [SNR_M2 ; snr_2 ] ; 

end 

ATR= ( 0 : .005:1) # ; 

res= [SNR_AR, SNR_M1 , SNR_M2 ] ; 

figure (16) 

plot (ATR, res, ATR, 3*ones (size (ATR) )  ,  'k:  ' )  ; 
axis([0  1  -20  max ( SNR_M1 ) ] ) ; 
xlabel ( 'AMPLITUDE  OF  ELVIS ' ) ; 
ylabel ( ' SNR  (dB)'); 

legend ( [ 'Whitening  Filter '],[ 'Method  1'], ['Method  2']); 
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B.  METHOD  1 


1.  Method  1  Filter  Process 

%Name :  J ohn  D .  Stevens 
%Date:  1  March  2000 

%Filename:  wavepack.m 

%description:  this  program  will  first  whiten  a  colored  noise 
signal 

%y  and  then  decompose  it  into  a  Haar  wavelet  decomposition  packet 
% framework  the  channels  will  be  white  noise  processes 
% 

function  [ STAGE 1 , STAGE 2 , STAGE3 , STAGE4 ] =wavepack (y ) 

if  size (y , 2 ) >size (y ,  1 )  ,  y =y ' ; end %make  the  signal  a  column  vector 

%perform  ARMA  process  to  get  a  whitened  signal 

[b, a] =cristi_ar (y) ; 

e0=filter (a,b,y) 

^filter  the  signal  through  the  Haar  wavelet  transform  packet: 

M=  [1/ sqrt (2) ,  1/sqrt (2 ) ; 1/1/sqrt (2 ) ,  -l/sqrt(2)]; 
temp=M* reshape (eO, 2 , length (eO) /2) ; 

Al=temp (1, : ) ; 

Dl=temp (2, : ) ; 

temp=M* reshape (Al, 2, length (Al) 12) ; 

AA2=temp (1, : ) ; 

DA2  =temp ( 2 ,  : )  ; 

temp=M*reshape {AA2 , 2 , length (AA2 ) /2 ) ; 

AAA3=temp (1,  :  )  ; 

DAA3=temp(2, : ) ; 

temp=M* reshape (DA2 , 2 , length (DA2 ) /2 ) ; 

ADA3  =temp ( 1 , : ) ; 

DDA3=temp (2 , : ) ; 

temp=M* reshape (Dl, 2 , length (Dl) /2); 

AD2=temp (1, : ) ; 

DD2=temp (2 , : ) ; 

temp=M*reshape (AD2 , 2 , length ( AD2 ) /2) ; 

AAD3  =temp ( 1 , : ) ; 

DAD3=temp ( 2 , : ) ; 

temp=M*reshape (DD2 , 2 , length (DD2 ) /2 ) ; 

ADD3=temp(l, : ) ; 

DDD3=temp (2 , : ) ; 

STAGEl=eO ; 

STAGE2  = [ Al ; Dl ]  ; 

STAGE3 = [AA2 ; DA2 ; AD2 ; DD2 ] ; 

STAGE4=  [AAA3  ;  DAA3  ;  ADA3  ;  DDA3  ;  AAD3  ;  DAD 3  ;  ADD3  ;  DDD3  ]  ; 
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%Name:  John  D.  Stevens 

%Date:  1  March  2000 

%Filename :  wavepackz .m 

%description:  this  program  will  first  whiten  a  colored  noise 
signal 

%y  and  then  decompose  ’it  into  a  Haar  wavelet  decomposition  packet 
% framework .  The  result  will  include  the  display  filtering 
process, 

% smoother 


function  [ STAGEl , STAGE 2 , STAGE3 , STAGE4 ] =wavepack (y ) 

if  size (y, 2 ) >size (y,  1)  ,  y =y 7 ;  end %make  the  signal  a  column  vector 

%perform  ARMA  process  to  get  a  whitened  signal 

[b,  a] =cristi_ar (y) ; 

e0=filter (a,b,y) ' ; 

%f liter  the  signal  through  the  Haar  wavelet  transform  packet: 

M= [ 1/sqrt (2 )  ,  1/sgrt (2 ) ; 1/1/sqrt (2 ) ,  -l/sqrt(2)]; 
temp=M* reshape (eO , 2 , length (eO) 12) ; 

Al=temp (1, : ) ; 

Dl=temp (2 , : ) ; 

temp=M*reshape(Al,2,  length (Al)  /2); 

AA2=temp(l, : ) ; 

DA2=temp(2, : ) ; 

temp=M*reshape (AA2 , 2 , length (AA2 ) / 2 ) ; 

AAA3=temp (1, : ) ; 

DAA3=temp (2 , : ) ; 

temp=M* re shape (DA2 , 2 , length (DA2) /2) ; 

ADA3=temp ( 1 , : ) ; 

DDA3=temp (2 , : ) ; 

temp=M*reshape(Dl/2/ length (Dl) /2) ; 

AD2=temp(l, : ) ; 

DD2=temp(2 , ; 

temp=M*reshape (AD 2 , 2 , length (AD2 ) /2 ) ; 

AAD3=temp (1, : ) ; 

DAD3=temp  (2  ,  : )  ; 

temp=M*reshape (DD2 , 2  # length (DD2 ) /2 ) ; 

ADD3=temp ( 1 , : ) ; 

DDD3=temp(2/ : ) ; 

STAGEl = smoother (eO) ; 

STAGE2=smoother ( [A1;D1] ) ; 

STAGE3= smoother ( [AA2 ; DA2 ; AD2 ; DD2 ] ) ; 

STAGE4= smoother  (  [AAA3  ;  DAA3  ;  ADA 3  ;  DDA3  ;  AAD3  ;  DAD 3  ;  ADD 3  ;  DDD3  ]  )  ; 
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2. 


Narrowband  Transient  SNR  Test 


%Kame :  John  D.  Stevens 

%Date :  11  Jan  2000 

%  F i 1 ename :  my_SNRl . m 

%description:  find  the  signal  to  noise  ratio  of  the  all  the 

channels 

%xn  the  Haar  Wavelet  Packet 
% 

%  SNR=my_SNRl ( STAGE1 , STAGE2 , STAGE3 , STAGE4 ) ; 

% 

function  SNR=my_SNRl ( STAGE1 , STAGE2 , STAGE3 , STAGE4 ) ; 

%the  assumption  is  that  the  transient  signal  is  .1  seconds  long 
NS_T=floor( .1*11025) ; 

Slp=l; 

Ns_tl=NS_T; 

Ns_t2=floor(NS_T/2); 

Ns_t4=floor(NS_T/4); 

Ns_t8=f loor (NS_T/8) ; 

%zfar  is  a  column  vector  while  all  the  other  inputs  are  row 
vectors 

%indexes  due  to  downsampling: 
indl=10000 ; 

indl_fin=indl  +  floor (Ns_tl*slp)  ; 

ind2=f loor ( indl/2 ) ; 

ind2_f in=ind2  +  floor (Ns_t2  *  s lp )  ; 

ind4=floor (indl/4) ; 

ind4_f in=ind4  +  f loor (Ns_t4*slp)  ; 

ind8=floor(indl/8) ; 

ind8_fin=ind8+f loor (Ns_t8*slp)  ; 

S=STAGE1 ; 

A1=STAGE2 (1, : ) ; 

D1=STAGE2 (2, : ) ; 

AA2=STAGE3 (1, : ) ; 

DA2=STAGE3 (2, : ) ; 

AD2=STAGE3 (3, : ) ; 

DD2=STAGE3 (4, :); 

AAA3=STAGE4 (1, : ) ; 

DAA3  =STAGE4 (2,  : ) ; 

ADA3 =STAGE4 (3,  : ) ; 

DDA3=STAGE4 (4, : ) ; 

AAD3  =STAGE4  ( 5  ,  :  )  ; 

DAD3=STAGE4 (6, : ) ; 

ADD3=STAGE4 (7, : ) ; 

DDD3=STAGE4 (8, : ) ; 
k=l .  5 ; 

nS=var ( [S (1 : (indl-1 )), S (floor (k*indl_f in) :length(S) )]); 
nAl=var( [Al(l: (ind2-l) ) , A1 ( floor (k*ind2_f in) : length (Al) )]); 
nDl=var( [Dl(l: (ind2-l) ) ,D1 (floor (k*ind2_f in) :length(Dl) ) ] ) ; 
nAA2  =var ( [AA2 ( 1 : (ind4-l) ) ,AA2 (floor (k*ind4_f in) : length (AA2) ) ] ) ; 
nDA2  =var ( [DA2 (1: (ind4-l) ) ,DA2 ( floor (k*ind4_f in) : length (DA2) ) ] ) ; 
nAD2  =var ( [AD2 (1: (ind4-l) ) ,AD2 ( floor (k*ind4_f in) :length(AD2) ) ] ) ; 
nDD2=var( [DD2 (1: (ind4-l) ) ,DD2 ( floor (k*ind4_f in) :length(DD2) ) ] ) ; 
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nAAA3=var  (  [AAA3  (1 :  (ind8- 

1)  )  ,AAA3  ( floor  (k*ind8__f in)  :  length  (AAA3 )  )  ]  )  ; 
nDAA3=var ( [DAA3 (1: (ind8- 

1)  )  , DAA3 (floor (k*ind8_f in) : length (DAA3 ) ) ] ) ; 
nADA3=var ( [ADA3 (1 : (ind8- 

1) ) ,ADA3 (floor (k*ind8_f in) : length (ADA3 ) ) ] ) ; 
nDDA3  =var { [ DDA3 ( 1 : ( ind8 - 

1) ) ,DDA3 (floor (k*ind8_f in) : length (DDA3 ) ) ] ) ; 
nAAD3  =var  (  [AAD3  ( 1 :  ( ind8  - 

1)  )  , AAD3 (floor (k*ind8_f in) : length (AAD3 ) ) ] ) ; 
n_DAD3  =var  (  [  DAD 3  ( 1 :  ( ind8  - 

1)  ) ,DAD3 (floor (k*ind8_f in) : length (DAD3 ) ) ] ) ; 
nADD3  =  var  (  [  ADD3  ( 1 :  ( ind8  - 

1)  ) ,ADD3 ( floor  (k*ind8__f  in) : length (ADD3 ) ) ] ) ; 
nDDD3 =var ( [ DDD3 ( 1 : ( ind8 - 

1) )  ,  DDD3 (floor  (k*ind8__f in)  : length (DDD3 ) ) ] ) ; 

%FIND  THE  SIGNAL  POWER: 
sS=var (S (indl : indl_f in) ) ; 
sAl=var  ( A1  ( ind2  :  ind2_f in)  )  ; 
sDl=var (Dl (ind2 : ind2_f in) ) ; 
sAA2=var  (AA2  (ind4  :  ind4_f  in)  )  ; 
sDA2  =var  ( DA2  (ind4  :  ind4_f  in)  )  ; 
sAD2  =var  ( AD2  ( ind4  :  ind4_f  in )  )  ; 
sDD2 =var ( DD2 ( ind4 : ind4_f in)); 
sAAA3=var  (AAA3  ( ind8  :  ind8_f  in)  )  ; 
sDAA3=var  (DAA3  (ind8  :  ind8_fin)  )  ; 
sADA3=var  (ADA3  ( ind8  :  ind8_f  in)  )  ; 

SDDA3  =var  (DDA3  ( ind8  :  ind8_f  in)  )  ; 
sAAD3=var  (AAD3  ( ind8  :  ind8_f  in)  )  ; 
sDAD3=var  (DAD 3  (ind8  :  ind8_fin)  )  ; 
s ADD 3=  var  (ADD3  (ind8  :  ind8_f  in)  )  ; 
sDDD3=var (DDD3 (ind8 : ind8_f in) ) ; 

%get  the  raw  ratio: 
snr_S=sS/nS; 
snr_Al=sAl/nAl  ; 
snr_Dl=sDl/nDl ; 
snr_AA2=sAA2  /nAA2  ; 
snr_DA2  =  s  DA2 / nDA2 ; 
snr_AD2=sAD2 /nAD2 ; 
snr_DD2=sDD2 /nDD2 ; 
snr_AAA3 = s AAA3  /nAAA3  ; 
snr_DAA3 =sDAA3  /nDAA3  ; 
snr_ADA3 = s ADA3  / nADA3  ; 
snr_DDA3  =  sDDA3  /nDDA3  ; 
snr_AAD3=sAAD3/nAAD3  ; 
snr_DAD3=sDAD3 /nDAD3 ; 
snr__ADD3  =  s ADD3  /nADD3  ; 
snr_DDD3=sDDD3 /nDDD3 ; 

c=  [ snr_S ;  snr_Al ;  snr_Dl ;  snr_AA2  ;  snr_DA2  ;  snr_AD2  ;  snr_DD2  ;  .  .  . 
snr_AAA3  ;  snr_DAA3  ;  snr_ADA3  ;  snr_DDA3  ;  snr_AAD3  ;  snr_DAD3 ;  snr_ADD3  ;  sn 
r_DDD3 ] ; 

x=c;x=x-ones(size(x) ) ; 

I=f ind (x< .01); 
x ( I )  =  . 01 ; SNR=10*logl0  (x)  ; 

SNR=  [ SNR  ( 1 )  ; max  (SNR  (2  : 15)  )  ]  ; 
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C.  METHOD  2 


1.  Method  2  Filter  Process 

%Name:  John  D.  Stevens 

%Date:  11  Jan  2000 

%Filename :  method_02_f ilt .m 

%description :  this  function  will  run  a  signal  through  the  first 

algorithm 

%that  was  developed  by  Dr.  Roberto  Cristi 
% 

%  [zl,zlt,z2,z2t,z3,z3t] ==method_02_f ilt (y) 

% 

function  [zl, zlt, z2 , z2t , z3 , z3t ] =method_02_f ilt (y) 

if  size  (y,  1)  >size  (y,  2)  ,  y=y' ;  end  %  y  row  vector 
k=l/sqrt (2) ; 

%k= . 5 ; 

[b, a] -cristi_ar (y ' ) ; 

[bl , al , bit , alt] =model2k (b,  a,  k) ; 

[b2 , a2 , b2t , a2t ] =model2k (bl , al , k) ; 

[b3 , a3 , b3t ,  a3t ] =model2k (b2 , a2 , k) ; 

M=k* [1 , 1; 1, -1] ; 

yl=M* reshape (y , 2 , length (y) /2 ) ; 
y2=M*reshape (yl (1, :)f2, length (yl) /2) ; 
y3=M*reshape (y2 (1, : ) , 2 , length (y2 ) /2) ; 

%e0=filter (a,b,y) ; 

el=f ilter (al , bl ,yl (1,:)); 

elt=f ilter (alt , bit ,yl (2 , : ) ) ;  %row  vector 

e2=filter (a2,b2,y2 (1, : ) ) ; 

e2t=f ilter (a2t , b2t ,y2 (2 , :)); 

e3=f ilter (a3 ,b3 ,y3 (1,  : )  )  ; 

e3t=f ilter (a3t,b3t,y3 (2,  : ) ) ; 

%z0=smoother (eO ) ; 
zl=smoother (el) ; 
zl=zl/max(zl) ; 
zlt=smoother (elt) ; 
zlt=zlt/max(zlt) ; 
z2=smoother (e2) ; 
z2=z2/max(z2) ; 
z2t=smoother (e2t ) ; 
z2t=z2t/max(z2t) ; 
z3 = smoother (e3 ) ; 
z3=z3/max(z3 ) ; 
z3t=smoother (e3t ) ; 
z3t=z3t/max(z3t); 
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2.  Multiresolution  Innovation  Filter  Bank  Routine 

%Name:  John  D.  Stevens 

%Date :  23  Feb  2000 

%Filename:  model2k.m 

%description :  This  program  applies  the  multiresolution 

innovation  process 
% 

function  [b2 , a2 , b2t , a2 t] =model2k (b, a, k) 

%  [b2,a2,b2t,a2t]=model2  (b,a) 

%  "t"  stands  for  "tilde" ,  which  is  the  high  pass  filter 
[Ce, Co, Fz] =polyd(b,a)  ; 

if  length(Ce) <length(Co) ,  num=k*Co+k* [Ce, zeros (1, length(Co) - 
length (Ce) )] ;  end 

if  length (Ce) >length(Co) ,  num=k*Ce+k* [Co, zeros (1, length (Ce) - 
length (Co) )] ;  end 

if  length (Ce) ==length (Co) ,  num=k*Ce+k*Co;  end 
den=Fz ; 

if  length (num) >length(Fz) ,  den= [Fz, zeros (1, length (num)  - 
length (Fz) )] ;  end 

[F1,G1,C1,D1] =tf2ss (num,  den) ; 

if  length(Co) +l<length (Ce)  ,  num=k*Ce+k* [0,Co, zeros (1, length(Ce) - 
length (Co) -1) ] ;  end 
if  length (Co) +l>length (Ce) , 

num=k* [ 0 , Co] +k* [Ce, zeros (1 , length (Co) +1 -length (Ce) ) ] ;  end 
if  length  (Co)  +l==length (Ce)  ,  num=k*Ce+k*-[  0  ,  Co]  ;  end 
den=Fz ; 

if  length (num) >length(Fz) ,  den= [den, zeros (1, length  (num)  - 
length (den) ) ] ;  end 

if  length (num) clength(Fz) ,  num= [num, zeros (1 , length (den) - 
length (num) ) ] ;  end 

[F2 , G2 , C2 , D2 ] =tf 2ss (num, den) ; 

%  overall  state  space 
if  length ( F2 ) ==length( FI) +1, 

Fl= [ [FI; zeros (1, length (FI) -1) ,1] , zeros (length (FI) +1,1)]; 

Gl= [G1 ; 0 ] ; 

Cl= [Cl , 0 ] ; 

end 


A=F1 ' ; 

B=[C2' ,C1'  ]  ; 
C=G2 ' ; 
D=[D2,D1] ; 

Q=B*B ' ; 

R=D*D ' ; 

S=B*D ' ; 


[K,  P]  =  mydlqe  (A,  eye  (size  (A)  )  , C,  Q, R,  S)  ; 
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[b2,a2] =ss2tf (A,K,C,  1) ; 
sigma=sqrt(C*P*C'+R) ; 
b2 =b2 /sigma ; 

%  "tilde"  -  the  high  pass  filtering 

if  length (Ce) < length (Co) ,  num=k*Co-k* [Ce, zeros (1, length(Co) - 
length(Ce) ) ] ;  end 

if  length (Ce)>length (Co) ,  num=k*Ce-k* [Co, zeros (1, length (Ce) - 
length (Co) ) ] ;  end 

if  length (Ce) ==length (Co) ,  num=k*Ce-k*Co;  end 
den=Fz; 

if  length (num) >length(Fz) ,  den= [Fz, zeros (1, length(num) - 
length(Fz) ) ] ;  end 

[F1,G1,C1,D1] =tf2ss (num,  den); 

if  length  (Co) +l<length  (Ce)  ,  num=k*Ce-k*[0,Co,zeros(l,length(Ce)- 
length(Co) -1) ] ;  end 

if  length (Co) +l>length (Ce) ,  num=k* [0,Co] - 
k* [Ce, zeros (1, length (Co) +l-length(Ce) ) ] ;  end 
if  length (Co) +l==length(Ce) ,  num=k*Ce-k* [0,Co] ;  end 
den=Fz ; 

if  length (num) >length (Fz) ,  den= [den, zeros (1, length (num) - 
length(den) ) ] ;  end 

if  length (num) <length (Fz) ,  num= [num, zeros (1 , length (den) - 
length(num) ) ] ;  end 

[F2 , G2 , C2 , D2 ] =tf2ss (num, den) ; 

%  overall  state  space 
if  length (F2 ) ==length (FI) +1, 

Fl= [ [FI; zeros (1, length (FI) -1 ) , 1] , zeros (length (FI) +1,1)]; 

Gl= [G1 ; 0] ; 

Cl= [Cl, 0] ; 

end 


A=F1 ' ; 

B= [C2 ' , Cl ' ] ; 

C=G2 ' ; 

D= [D2,D1] ; 

Q=B*B' ; 

R=D*D' ; 

S=B*D' ; 

[K,  P]  =  mydlqe  (A,  eye  ( size  (A)  )  ,  C,  Q,  R,  S)  ; 

[b2t,a2t] =ss2tf (A, K,C, 1) ; 
sigma=sqrt (C*P*C'+R) ; 
b2 t=b2t/ sigma; 
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3, 


Narrowband  Transient  SNR  Test 


%Name :  John  D.  Stevens 

%Date :  11  Jan  2000 

%Filename :  my_SNR2 ,m 

%description:  this  function  will  determine  the  SNR  of  various 

%transients  in  white  noise 
% 

%  SNR=my_SNR2 (y) 

% 


function  SNR=my_SNR2 (y) 

if  size (y, 1) >size (y,  2) ,  y=y';  end  %  y  row  vector 
kk=l/sqrt (2)  ; 

%kk^ .  5 ; 

[b,a] =cristi_ar (y' ) ; 

[bl,al,blt,alt]  =model2k(b,a,kk)  ; 

[b2/a2/b2t/a2t]  =model2k  (bl ,  al ,  kk)  ; 

[b3 ,  a3  ,  b3 1 ,  a3 1  ]  =model2k  (b2 ,  a2 ,  kk)  ; 

M=kk* [1,1;!, -1] ; 
yl =M* re shape  (y  ,2 ,  length  (y )  12)  ; 
y2=M* reshape (yl (1, : ) , 2 , length (yl) /2) ; 
y3=M*reshape  (y2 (1, : ) ,2, length (y2) 12)  ; 

%e0=filter (a,b,y) ; 

el=f ilter (al,bl,yl (1,:)); 

elt=f ilter (alt , bit ,yl (2 , : ) )  ;  %row  vector 

e2=f ilter (a2 , b2 ,y2 ( 1 , :)); 

e2t=filter (a2t,b2t,y2 (2,  : ) ) ; 

e3=f ilter (a3,b3,y3 ( 1 , : ) ) ; 

e3t=filter (a3t,b3t,y3 (2,  : )  )  ; 


%find  SNR: 

NS_T=floor( .1*11025) ; 
slp=l; 

Ns_tl=NS_T; 

Ns__t2  =  f  loor  (NS_T/2  )  ; 

Ns_t4=f loor (NS_T/ 4 ) ; 

Ns_t8=f loor  (NS_T/8 )  ; 

%zfar  is  a  column  vector  while  all  the  other  inputs  are  row 
vectors 

%indexes  due  to  downsampling: 
indl=10000; 

indl_f  in=indl+f  loor  (Ns__tl*slp)  ; 

ind2=floor (indl/2) ; 

ind2_f in=ind2  +  f loor (Ns_t2*slp)  ; 

ind4=f loor (indl/4) ; 

ind4_f in=ind4+f loor (Ns_t4*slp)  ; 
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ind8=floor (indl/8)  ; 

ind8_f in=ind8+f loor (Ns_t8*slp) ; 

zf 0=el; 
zf Ot=elt; 
zf l=e2 ; 
zf lt=e2t ; 
z  f 2=e3 ; 
zf 2t=e3t ; 

k=l . 5 ; 

%find  the  variance  of  the  noise 

mnO=var  (  [zfO(l:  (ind2-l)  )  ,  zfO  (floor <k*ind2_f in)  rlength(zfO)  )  ]  )  ; 
mnOt=var( [zfOt(l: (ind2- 

1) ) , zfOt (floor (k*ind2_f in) : length (zfOt) ) ] ) ; 

mnl=var ( [zf 1 (1 : (ind4-l) ) , zfl (floor (k*ind4_f in) : length (zfl) ) ] )  ; 
mnlt=var ( [zflt (1 : (ind4- 

1) ) , zflt (floor  (k*ind4__f  in) : length (zflt ) ) ] ) ; 

mn2=var ( [zf2 (1: (ind8-l) ) , zf2 (floor (k*ind8_f in) : length (zf 2) ) ] ) ; 
mn2t=var  (  [zf2t  (1:  (ind8- 

1) ) / zf2t ( floor (k*ind8_f in) : length (zf2t) ) ] ) ; 

%find  the  variance  of  the  spike  in  the  region: 

mpO=var (zf 0 (ind2 : ind2_f in) ) ; 

mpOt=var (zfOt (ind2 : ind2_fin) ) ; 

mpl=var  (zfl  (ind4 :  ind4__f  in)  )  ; 

mplt=var (zflt (ind4 : ind4_fin) ) ; 

mp2=var ( zf 2 (ind8 : ind8_f in) ) ; 

mp2t=var (zf2t (ind8 : ind8_fin) ) ; 

%get  the  raw  ratio: 
snr_0=mp0/mn0 ; 
snr_0t=mp0t/mn0t ; 
snr_l=mpl/mnl ; 
snr_lt=mplt/mnlt ; 
snr_2=mp2/mn2 ; 
snr_2 1  =mp2 1 / mn2 1 ; 

%SNR=  [  snr_ar ,  snr_0 ,  snr_0t ,  snr_l ,  snr_lt ,  snr_2  ,  snr_2t]  ; 

SNR=  [  snr_0 ,  snr_0 1 ,  snr_l ,  snr_l  t ,  snr_2 ,  snr_2 1  ]  ; 

%  check  for  ill  conditions 
X=SNR; 

x=x-ones (size (x) ) ; 

I=find(x< .01)  ; 
x ( I ) = . 01 ; 

SNR=10*logl0 (x) ; 

SNR=max  ( SNR )  ; 
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